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Hopkins  University  Applied  Physics  Laboratory,  which  greatly 
improved  the  accuracy  of  the  confidence  intervals  for  CEP. 
The  purpose  of  this  thesis  is  to  develop  and  test  procedures 
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I .   INTRODUCTION 

A.   BACKGROUND 

In  1966  W.R.  Blischke  and  A.H.  Halpin  [Ref.  1]  first 
described  methods  for  the  computation  of  confidence  intervals 
for  circular  error  probability  (CEP)  based  on  first  order 
variance  estimates.  When  these  procedures  were  implemented  at 
the  Johns  Hopkins  University  Applied  Physics  Laboratory  (JHU- 
APL)  it  was  found  that  under  certain  conditions  the  resulting 
confidence  intervals  for  CEP  were  smaller  than  expected.  As 
a  result  R.C.  Ferguson,  K.V.  Kitzman  and  P.B.  Jackson  [Ref.  2] 
have  extended  the  Blischke-Halpin  procedures  to  create  a 
second  order  variance  estimate.  Testing  conducted  by  Kitzman 
[Ref.  3]  has  confirmed  that  this  method  produces  a  more 
accurate  estimate  for  the  variance  of  CEP. 

Our  goal  will  be  to  extend  this  procedure  to  the  3- 
dimensional  case  to  obtain  a  second  order  estimate  for 
variance  of  spherical  error  probability  (SEP) . 

B.   METHODOLOGY 

We  begin  by  assuming  that  missile  detonations  are 
distributed  as  trivariate  Gaussian.  In  the  past,  analysis  of 
ballistic  missile  test  firing  data  has  failed  to  disprove  this 
assumption.   If  the  mean  and  variance  are  given  by: 


v>  = 


1*1 

^31 


E  = 


°11 

ai2 

°13 

°21 

°22 

°23 

°31 

°32 

°33. 

Then    the    probability    of    a    detonation    occurring   within    a 
distance  R  of   the   origin   is   given  by: 


d     (v2ti)  q         v  J 


where, 


x  = 


X, 


g2  =  det[2], 


and    £>  =  {(x1,x2,x3)|x12+x22+X32  £  i?2} 


.     .  .  dp 

By  the  implicit  function  theorem,     -5-   >  0   implies  that  the 


equation  P{R;\i,H)  =  CONSTANT  defines  a  dif f erentiable  function 


2?(li,E)  such  that  P(R(\x,  2);  \i ,  S)  =  CONSTANT.   If  the  CONSTANT  is 
set  to  — ,  then  i?(n,E)  is  the  SEP  function.   Let  Xx,Xz,...,Xn   be 


a  random  sample  of  independent  and  identically  distributed 
random  vectors  with  distribution  N(\i,2,)  .       Then  the  maximum 

1  A 
likelihood   estimators   of  \x      and   S   are   p  =  —  ^Xd      and 


2  =  —  ^2(Xi-ii)T(Xi-ii).      When  jl  and  2  replace  ja  and  S  in  the 

ni=l 


equation  £1(.R;|i,,E)  =  —  then  £((1,2)  gives  SEP,    an  estimator  of 


SEP. 

The  second-order  approximation  for  the  variance  of  SEP   is 
given  by: 

°1*p  =  (D^uSEP)TP(DiiJ,uSEP) 


+  ^(d^sep)  T(P®P)(d^uSEP) 

\[{d^sep)te]2  (1) 


+    — 


where  Su  =  (allf  o12,  o13,  o22,  a23,  o33)r/   D    ZuSEP  is  the  derivative  of 
SEP  with  respect  to   n.,Eu  viewed  as  a  column  vector,    D^USEP  is 


the  usual  second  derivative  matrix  (or  Hessian)  of  SEP  with 


respect  to  H,2U,  P  = 


\P\xZu     P2U ) 


is  the  variance-covariance 


matrix  ji,  2U,  ®     is  the  tensor  product  operation  and  the 

underline  notation  is  used  to  represent  the  column  vector 
formed  by  stringing  out  the  elements  of  the  matrix  by  rows. 


For  example,  if  A   = 


1  2 

3  4 


then  A  =  [1   2  3  4]T.   Note  that  (I  is 


1  S 

distributed  N(a,—  E)  so  P..  =   — .   It  can  also  be  seen  that/ill 

n  ■*    n 


has  a  Wishart  distribution  with  parameter   2  ,   so  that 
Ps  =  —  (202)  .   The  complete  derivation  of  equation  (1)  ,  based 


on  Taylor  Series  expansion,  is  given  in  appendix  A. 


Since  \x   and  2  are  independent  [Ref.  4:p.  102]  it  follows 
that  P  /Su  =  0  and  equation  (1)  can  be  written  as: 


„2 

°sSp 


'  dSEPYp(dSEP\ 


i^J 


J  dSEPyp    I  dSEP\ 


au2  J  {  "     >} 


(d*_SEP\J&SEP 


(P2U®P2U) 


'  &SEP 

v  as"2 


(2) 


Derivation  of  the  first  order  terms  in  equation  (2)  has 
been  accomplished  by  S.D.  Hill  at  JHU-APL  [Ref.  5].  In  order 
to  solve  for  the  second  order  terms  it  will  be  necessary  to 


find  expressions  for 


&SEP         PSEP 


d\i: 


asu; 


and 


&SEP 
d\id?,u' 


i.e.  it  is 


necessary  to  find  D2  j,uSEP ,  the  second  derivative  matrix  of  SEP 

with  respect  to  its  arguments.   Chapter  II  will  describe  the 
derivation  of  the  second  order  variance  estimate  and  Chapter 
III  will  discuss  implementation  and  testing  of  the  algorithm. 
Conclusions  will  be  presented  in  Chapter  IV. 


II.   DERIVATION  OF  THE  SECOND  ORDER  ESTIMATE 

A.   FIRST  ORDER  TERMS 

We  begin  by  describing  the  derivation  of  the  first  order  terms 
given  in  eguation  (2) .   If  we  define: 


where 


1-1  . 


a11 

a12 

o13' 

a21 

a22 

a23 

a31 

o32 

a33. 

o12  =  o21 
a13  =  a31 
o23  =  a32 


(x-ia)^-1^-^  (x1-^1)2o11+2(x1-^1)(x2-^2)o21+2(x1-^1)(x3-^3)a31 
+(x2-ti2)2a22+2(x2-n2)(x3-^3)o23+(X3-n3)2a33 


and 


p(p,e,<f>;  H,  Su)  =  i^psinecosct),  psin6sin<J>,  pcos8;  \i,  Su) 


Then,  with  a  shift  to  spherical  coordinates,  we  see  that: 


P(i?;H,Eu)  =  ^—  f  *('(   p2sin6g<p,e,<t);^/Eu)c?pded(J) 

l.rT^VrrJ0     JO  JO 


(/2T)Jg 


(3) 


We  wish  to  solve  equation  (3)  for  each  of  the  partials 

dSEP 


axi 


where  X   =  (\ix,  \i2,  \i3,  o11,  o12,  o13,  a22,  o23,  o33) .     From   the 


previous  section,  i?(n,E)  =  — ■  is  the  SEP  function,  so  we  need 


dR 
to  solve  for  the  partials  -^r— . 


Taking  partials  of  equation  (3)  we  see  that, 


and 


P'(R)   =  — f  *  f"sinBg(R,e,$)ddd<\>. 

(v/2¥)V°  Jo 

We  can  solve  the  right  hand  side  of  equation  (4)  for 
i=l,2,3  by  using  the  relationships 

_dg_=__dg_  J>9_=-_d9_  dg  __  dg 

d\xx       dx1  d\iz       dx2  d\i3       dx3 


so  we   see   that 


-/// 


-&dxtdx2dx3 


■III 


ex,      L     *     J 


If  we  apply  Green's  Theorem  to  this  equation  and  then  make 
a   shift  to   spherical   coordinates  we  get: 

-fff-&dx1dx2dx3  =  R2f2*j*cos<bsin26g(R,Qf<b)dBd<b 

D  1 

-fff-^-dXidXzdXi  =  R2f2rcj*sin$8in26g(R,B,$)dQdto 


- j j f  -& dxxdx2dx3  =  R2[2"  r co s8sin6 g(R,  8,  <t>)d0d<J> 

D  3 

=    *L-f2*[Ksin2Qg(RlQl<b)dQd<b 
2  Jo    Jo 


To  solve  the  right  hand  side  of  equation   (4)    for  i=4 ,  ...  ,9 
we   first  note  that: 


&9       _  , 

do,.    ' 

i=j 

d\Lxd\L2 

dg 
da..    ' 

i*j 

so  we  get 


dg 

=  \ 


do±J 


2  dxj  dpi   '  J 


which   is   used   to   evaluate   the    integral   of    -rrf-    as    follows 


- f  f  f -$?- dx1dx2dx2  = — f     f  [a11(psin0cos4>-jjL1)+o12(psin0sin4)-^2) 

D  4 

+  a13(pcos6-^3)]cos<t)sin20g(^/  8,  4>)d0d(J) 

-  I  f  f  -^- dx1dx2dx2  =  R2  f  'tr't[a11(psin8cos<|)-^1)+o12(psinesin<|)-n2) 

D  5 

+a13(pcos0-^3)]sin<|)sin20g^i?,  6,  <J))d0d(J) 

-  f  j  I -Jr-  dx1dx2dx2  =  —\      j    [a11(psin0cos(j)-ji1)+o12(psin0sin(J)-ji2) 

D  6  ° 

+a13(pcos0-^3)]sin20gt^,  6/  4>)d0d<t> 

-  \  \  \ -^- dx1dx2dx2  =  ——  j      j    [o21(psin0cos<t>-ji1)+o22(psin0sin<|)-jjL2) 

d         7 

+  a23(pcos0-|i3)]sin<t>sin20g(i?,  0,  (J>)d0d<t> 


-  j  (  j  -^f-  dxxdx2dx2  = — j      j    [a21(psin0cos(|)-^ji1)+o22(psin0sin(l)-n2) 

D  8  0         0 

+a23(pcos0-n3)]sin20g(i?,0,<J))d0d<t) 


-  f  f  f  _HjL  dxxdx2dx3  =  —  r2*rn[^31(psin0cos<|)-n1)+o32(psinesin<j)-ji2) 

D  9 

+o33(pcos9-n3)]sin26g(i?/  6,  <$)dQd<b 


If  we  define  the  functions: 


CCd   ,  k   ,(^;n,Su)  =  f   " [Kcosi(jQ)cosk(li\>)g{Rldl^)d6d^ 

' J  '     '  J  o  J  0 


CSiijikil(R;\i,Xu)   =  |  Y*cosiO'0)sinV*)9(tf,0,  *)<#<*!> 


SCi#i#Jtfl(R;|i#E,0  =  f  1tf%sini(jB)cosk(l$)g(R,ei$)dBdb 


ssi.j.k.jiRiP>Eu)  =  f  n (nsini(jQ)sin\li>)g{Rief^)ded^ 


Extending  the  work  of  Blischke  and  Halpin,  with  the  aid  of 
Green's  Theorem,  we  can  see  that: 


dp. 


SC2  ,1,1,1   "^2 ,1,1,1    "^^1 ,2,0,0 


Of,  qp  p  q(~, 

ol~l, 1,0,0   °*-i,  1,0,0   ^°°1, 1,0,0. 


(5) 


l(l<M.Su);|».IIu) 


10 


dSEP(\i ,  Su) 


2SC 


1,1,0,0 


dSC2.l.l.l 


d\i1 


dSS. 


2,1,1,1 


d\xx 

oSCi  ,2,0,0 


d\i1 


dSS. 


2,1,1,1 


d\l- 


dSC. 


1,2,0,0 


d\i2 

vSCi  ,2,0,0 


2d\x. 


(6) 


l(j<^.  xu);n,i:u). 


Applying  the  partials  given  in  eguation  (6)   to  the 
partials  - f  j  J  -5^- dx1dx2dxi   given  above  we  see  that: 


dSEP 


1    [R 


-^8(°11JJll+0l2^2+al3^3)} 


'     =     -^(t[2011^1+4°12^+01>3-^)] 


dax. 


w5  1  4 


Vlo(oll^l+«yl2^2+alV3)} 


dSEP 
da12 


1    f*r„i 


2wc 


-^(0Uh+012f»2+013h)} 


11 


dSEP   =      1    [R 
do22  2w5 


||[2o12^11+4o22^12  +  o2>3-^)] 
-v10(o12ji1+o22ji2+a23ji3)} 


dSEP  1    [R 


do23  2w5\ 


[a12^w1-w2)+o22(wz-w4)+a23^w1+w5)j 


-V6(o12^1+o22ji2+o2V3)} 


-^6(o13Ji1+o23n2+a33ji3)} 


do33  4w5 


where    fo{i?;n,Eu)    is   given  by: 


^,^,2U)  = 


wi  "  ^-W,  1,1,1 
w2   ::  CC1311 

W3     '"    *-''->  1 , 1 , 1 , 1 
Wi    ~    ^1,3,1,1 

V6     =    ^-"1,2,0,0 
^7     "    ^^-1,3,0,0 

w8  •-  SC2tX11 
w9   :=  SC3>121 


w. 


SS2l  1,1,1 


vii   =  SS3,i,i,: 

^12     "    ^^3 ,1,2, 


where  each  of  the  functions  is  evaluated  at  R(\i,  Eu),  \i,  Eu . 

This  completes  the  derivation  of  the  first  order  terms 
required  in  the  computation  of  the  variance  estimate. 


12 


B.   CONCEPT  FOR  DERIVATION  OF  THE  SECOND  ORDER  TERMS 

We  first  observe  that  equations  (5)  and  (6)  enable  one  to 

express  D*  ^USEP      as   a   composition   of   three   functions 

Tlt    T2   and  T3 ,    each  of  which  is  a  fairly  simple  function  of  its 

arguments. 

Explicitly, 


where, 


r>,Eu)  = 


SEP(\x,Zu) 


r>,Eu,i?)  = 


R 


13 


Vc 


We 


w. 


10 

-^(4o11v9+2a12Vll+a1>1-v2))-^(a1^1+o12ji2+o13ji3) 

^1 


4w, 


(2o11^11+4a12v12+a13(V3-w4))--^(a11^1+o12^2+a13n3) 


w^ 


R 


w, 


Awc 


(o1\w1-w2)+o12(w3-Wt)+o12(w1+w5))--^(o11\x1+o12\i2+o13\i3) 


2  Wc 


R 


V, 


v-(2o12w11+4o22Vl2+a2>3-^))-^(a12ji1+o22^2+o23^3) 


E 


w. 


v•(o12(Vl-^2)+o22(^3-^)+o23(w7+V5))-^-(a12^i1+o22^l2+a23^3) 


Wc 


8v, 


(o13(^1-v2)+o23(^3-^4)+a33(w7+v5))-— ^-(a13ji1+o23ji2+a33n3) 


4  We 


Since  D  ZuSEP(\i ,  V)  =  r3°r2or1(|i/ Su)  we  may  apply  the  chain 
rule   to   obtain  the   second  derivative   as   a  matrix  product: 

d\  ^sEPiii ,  s u)  =  r>r3(r2«>r> ,  s u))  -Dr^fo  ,  E u))  -dt^  ,  E ") .        ( 7 ) 

The  simplicity  of  the  functions  rx,  r2  and  r3  means  that 
the  derivatives  DTX,  DT2  and  DT3  are  easily  obtained  and  then 
the  rather  large  matrix  products  involved  in  DT3-DT2'DTX  can  be 
left  to  the  computer  leading  to  a  rather  painless  computation 

of  D2  yuSEP. 


14 


DTX     (the  total  derivative  of  the  function  T1)    is  the  10x9 
matrix  consisting  of  the  9x9  identity  matrix  (from  -31L 1\ 


with  the  tenth  row  consisting  of  the  partials  of  SEP  with 
respect  to  mu  and  sigma. 


DTX   = 


0 
0 
0 
0 
0 
0 
0 
0 

1 

dSEP     dSEP     dSEP     dSEP     dSEP     dSEP     dSEP     dSEP     dSEP 
d\ix        d\i2        d\x2       da^      da12      do12 


1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

da22      5o23  da33 


DT2    is  the  22x10  matrix  consisting  of  the  10x10  identity- 
matrix  (from  ^^  ' — '-)  ,    with  rows  11  to  22  consisting  of  the 


12x10  sub-matrix  of  partials  of  W  with  respect  to  [i,    Eu  and 
R. 
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1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

DT2   = 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

dwx 

dw1 

dw1 

dwx 

dwx 

dwx 

dw1 

dwx 

dw1 

dw1 

d\i1 

d\i2 

d\i3 

d°ll 

d°12 

da13 

da22 

da23 

da33 

dR 

dw12 
d\i1 

dw12 
dp2 

dw12 
d\i3 

dw12 
daxl 

3ai2 

dw12 
da13 

dwl2 

d<J22 

dw12 

3°23 

3V12 

do33 

dw12 

dR 

DT3    is  the  9x22  matrix  consisting  of  the  partials  of 
ui(i  =  l, ....  9)  with  respect  to  \i,    Eu,  R  and  W. 


DT3   = 


du,  du,  3u,  du. 

■^(1x3)   ^i(lxS)   -^(lxl>   ^(1X22, 


dus 

"ail 


(1X3) 


3E" 


(1x6) 


du9 
~8R 


(1x1) 


Bu9 


(1x12) 


Thus,  in  order  to  obtain  the  second  degree  variance 
estimate  given  in  equation  (2)  we  must  solve  for  the  values  of 
all  elements  of  DTX,   L>T2  and  Dr3 .    Once  these  values  are 

computed,  the  solution  is  found  by  the  matrix  multiplication 
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shown  in  equation  (7)  ,  yielding  the  9x9  matrix  of  second 
partials  of  SEP  with  respect  to  mu  and  sigma. 

The  values  of  the  DTX   elements  are  computed  from  equations 

(5)  and  (6)  .  The  determination  of  DT2   requires  the  evaluation 

of  -^— ,  and  -rr-   which  is  fully  described  in  appendix  B. 

d\x       asu      BR 

Finally,  to  solve  for  the  elements  of  DT3    it  is  necessary  to 

solve  for  each  of  the  partials  given  above.  The  evaluation  of 
these  partials  is  given  in  appendix  C. 
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III.   IMPLEMENTATION  AND  TESTING 

A.   PROGRAMMING 

The  calculation  of  o2sgp   the  second  order  approximation  for 

variance  of  SEP  (when  \i   and  S  are  known)  has  been  implemented 

in  the  FORTRAN  program  SEPCMP  which  is  provided  in  appendix  D. 

The  following  steps  are  followed  to  obtain  the  approximation: 

Input  the  distribution  mean  (\i)      and  variance- 
covariance  matrix  elements  (Su) . 
Compute  SEP  using  \x    and  S  . 

-  Compute  the  values  of  each  element  of  the  matrices 
DTI,  DT2  and  DT3  and  then  form  the  matrix  product 
comprising  the  second  derivative  of  SEP  with 

respect  to  \x  and  Eu.  The  first  derivative  of  SEP 
with  respect  to  \x  and  Su  is  then  given  by  row  10 
of  the  matrix  DTI. 

-  Calculate  the  variance-covariance  matrices  of  (1 

(P^)  and  of  2  (PE)  from  2. 

Form  the  SEP  variance  approximation  using  equation 
(2). 

Computation  of  SEP  given  \x    and  Su  is  based  on  the  PL-1 

program  RAP  provided  by  JHU-APL.   RAP  is  a  general  program 

which  provides  the  radius  of  coverage  for  any  probability  in 

the  interval  (0,1)  for  either  2  or  3  dimensions.    Those 

portions  of  the  program  applicable  to  this  problem,   ie. 

probability  =  0.5  and  3  dimensions,  have  been  used  in  the 

FORTRAN  subroutine  FINDSEP  to  obtain  SEP.   Derivation  of  the 

equation  used  to  compute  SEP  is  given  by  L.  S.  Simpkins  [Ref. 

6]  and  is  implemented  in  the  subroutine  FINDSEP  using  Gaussian 
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Quadrature  (via  the  IMSL  subroutine  QTWODQ)  to  evaluate  the 
required  double  integral. 

Gaussian  Quadrature  is  also  used  to  evaluate  each  of  the 
double  integrals  in  the  W,  CC,  CS,  SC  and  SS  functions,  which 
are  required  in  the  evaluation  of  the  matrices  DTI,  DT2  and 
DT3. 

B.   TESTING 

The  purpose  of  deriving  the  second  order  approximation  of 

the  variance  of  SEP  is  to  yield  more  accurate  confidence 
intervals  for  SEP.  The  approximation  of  these  confidence 
intervals  is  derived  from  the  fact  that  since  p.  and  H  are 
maximum  likelihood  estimators  of  \x  and  S  and  thus  SEP  is 
asymptotically  Normal  with  mean  SEP.  Thus,  for  example,  an 
approximate  95%  confidence  interval  for  SEP  is  given  by 

SEP  ±1.96  yJvar(SEP) . 

In  order  to  test  the  affect  of  using  the  second  order 
approximation  for  the  variance  of  SEP,  tests  of  confidence 
interval  coverage  were  conducted  as  follows: 

Sample  3  0  random  vectors  from  i\^p,,E)  and  compute 

|i  and  U  for  the  sample. 
Calculate  SEP   from  ji  and  2  . 

-  Calculate  o2sip   by  replacing  £  and  2  for  n  andE 
in   both   the   first   order   and   second   order 

approximations  for  the  variance  of  SEP. 

-  Compute   the   large   sample   approximation   95% 

confidence  intervals  for  SEP  using  SEP±1 .  9  6^6^. 
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Repeat  1000  times  to  calculate  the  percentage  of 

times  that  confidence  intervals  based  on  each 

approximation  cover  the  true  SEP. 

For  each  input  variance-covariance  matrix,   S , 

repeat  the  test  at  different  magnitudes  of  the 

mean. 

Trivariate  normal  random  vectors  are  obtained  using  the 
IMSL  subroutines  DCHFAC  and  DRNMVN.  DCHFAC  performs  Cholesky 
factorization  of  the  distribution  SIGMA  which  is  then  used  as 
an  input  to  DRNMVN  which  returns  the  desired  random  vector. 

The  estimates  of  P^    and  PE  are  found  by  replacing  \x    and 

2  with  |2  and  2  in  the  formulations  of  P^    and  Ps  given  in 

Chapter  I.   The  results  of  this  test  are  summarized  in  Table 
1. 

TABLE  1.   95  PCT  CONFIDENCE  INTERVAL  COVERAGE 


SIGMA 

MU 

FIRST  ORDER 
CI 

SECOND  ORDER 
CI 

1.0     0.8     0.9 
0.8     1.0     0.8 
0.9     0.8     1.0 

0 
0 
0 

94.6 

95.2 

1.0 
1.2 
0.8 

92.8 

93.3 

100    -160     270 

-160     400    -480 

270    -480     900 

0 
0 
0 

93.3 

94.2 

100 
250 
500 

95.6 

95.8 

10,000   3,000   2,000 
3,000  10,000   3,000 
2,000   3,000  10,000 

0 
0 
0 

93.3 

93.8 

8,000 
10,000 
12,000 

93.5 

94.0 
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These  results  indicate  that  there  is  little  accuracy  to  be 
gained  by  using  the  second  order  variance  approximation. 
There  is,  however,  a  situation  in  which  the  second  order 
approximation  becomes  important.  This  case  has  been  described 
by  K.  V.  Kitzman  [Ref.  3]  and  is  paraphrased  here. 

Weapons  system  targeting  can  be  measured  from  separate 
subsystems  (such  as  navigation,  guidance,  postboost,  etc.). 
Suppose  that  there  are  n  independent  subsystems  and  let 

Xilf...,Xim   ~  N(\x,11)    represent  the  errors  for  the  1th    subsystem 
i  =  1, . . . ,n,  then  the  impact  errors  for  the  weapons  system  are 

n  i  iH 

Yj   =  X)*ij  j=l,...,m  and  Ylt...,Ym   ~  N(\xY=n\il'LY=n'L) . 

2=1 


If  the  Y's  are  measured  directly  and  used  to  estimate \x 


Y 


2 

and  Sy   (to  get  the  system  SEP)   then  P     =  — -  =  —  2  and 

v'y        m         m 


psy  =  —  (sr®sr)  "  -^-(2®2).   If  on  the  other  hand  the  individual 


X's  are  available  from  the  subsystems  then  (iy  =  ^  X±    where 

i=l 


*i  =  ^txij   and  *r  =  iE(xiJ-NT(xij-^y)-        Then  p»r  =  ^ 


and 
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P2  =  —  (2®E).  The  significant  factor  is  that  the  estimate  of 


Ps  differs  from  the  first  estimate  by  a  factor  of  n  while  the 

estimate  of  PM   remains  unchanged.   Thus  as  the  number  of 

subsystems  becomes  larger  and  the  number  of  observations 
becomes  smaller  the  size  of  Ps  can  become  very  small  in 

relation  to  P^ .   As  the  mean  approaches  0  it  follows  that 

dSEP 


dH 


approaches  0  as  well  and  thus  the  contribution  of  P  to 


the  estimate  can  become  insignificant  even  though  it  is  large 
in  relation  to  the  PE  term. 

This  affect  has  been  examined  as  follows: 

-  Divide  \i   and  E  by  n,  the  number  of  iid  subsystems 
being  simulated. 

Draw  m  random  samples  from  the  distribution  given 
by  this  new  \x  and  E  for  each  of  the  n  subsystems. 
Calculate  p.  and  E  for  each  of  the  n  subsystems. 

-  Sum  the  n  estimates  to  get  a  total  (I  and  E  ,  the 
estimate  for  mean  and  variance  in  Y. 

Compute  confidence  intervals  for  SEP  using  both 
the  first  and  second  order  approximations  with 

Pu  =  -E  and  PE  =  —  (E®E). 


This  test  was  conducted  with  n  =  4  and  m  =  5.  The  results 
of  this  test  are  given  in  Table  2.  As  can  be  seen,  the  effect 
is  greatest  when  the  mean  is  0,  and  decreases  as  the  mean 
becomes  larger,  which  is  what  would  be  expected. 
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TABLE  2.   95  PCT  CI  COVERAGE,  MULTIPLE  SUBSYSTEMS 


SIGMA 

MU 

FIRST  ORDER 
CI 

SECOND  ORDER 
CI 

100    -160     270 

-160     400    -480 

270    -480     900 

0 
0 
0 

92.4 

97.6 

5 
10 
20 

91.4 

94.8 

25 

50 
80 

89.4 

92.2 

50 
150 
200 

92.6 

94.4 

To  test  the  sensitivity  of  the  second  order  variance 

approximation  procedure  to  the  assumption  that  the  detonations 

are  normally  distributed,  tests  were  conducted  in  which  random 

vectors  were  drawn  from  a  mixed  normal  distribution.  The  test 

is  conducted  in  the  following  manner: 

Prior  to  each  sample  obtain  random  variable  p  from 

£7(0,1)  distribution. 

If  p  is  less  than  .98  then  draw  a  sample  from 

N(\i,2,)    and  proceed. 

if  p  is  greater  than  .98  then  draw  a  sample  from 

N{\i, 10*1,)    and  proceed. 

Continue  until  3  0  samples  have  been  drawn,  then 

compute  all  values  and  proceed  as  above. 

This  procedure  simulates  a  distribution  with  the  same  mean 

but  larger  tails  than  the  normal  distribution.  The  confidence 

interval  coverages  are  then  compared  to  those  obtained  in  the 

base  case  from  the  same  input  distribution.   The  results  of 

this  test  are  summarized  in  Table  3. 
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TABLE  3. 

95  PCT 

CI  COVERAGE,  MIXED  DISTRIBUTION 

SIGMA 

MU 

FIRST  ORDER 
CI 

SECOND  ORDER 
CI 

1.0 
0.8 
0.9 

0.8 
1.0 
0.8 

0.9 
0.8 
1.0 

0 
0 
0 

86.6 

87.6 

1.0 
1.2 
0.8 

91.2 

91.8 

100 

-160 

270 

-160 

400 

-480 

270 

-480 

900 

0 
0 
0 

88.2 

88.4 

100 
250 
500 

91.5 

92.4 

10,000 
3,000 
2,000 

3,000 

10,000 

3,000 

2,000 

3,000 

10,000 

0 
0 
0 

79.3 

80.4 

8,000 
10,000 
12,000 

85.6 

86.2 
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IV.   CONCLUSIONS 

Examination  of  the  results  in  Table  1  shows  that  the 
second  order  approximation  procedure  provides  a  slight 
improvement  in  the  95%  confidence  interval  coverage.  The  real 
utility  of  this  method  is  brought  out  by  the  results  of  the 
multiple  subsystem  test  shown  in  Table  2 ,  where  confidence 
interval  coverage  improves  by  as  much  as  5%.  Since  this  case 
represents  the  actual  situation  in  missile  testing  (i.e.  very 
small  mean  and  independent  system  observations)  it  is 
reasonable  to  state  that  use  of  this  procedure  will  improve 
the  estimate  of  the  variance  in  SEP.  Furthermore  it  seems 
highly  unlikely  that  there  would  be  any  advantage  to  be  gained 
by  pursuing  a  higher  order  estimate  of  the  variance  of  SEP. 

The  mixed  distribution  test  shows  that  the  second  order 
approximation  is  only  slightly  more  robust  than  the  first 
order  approximation  if  the  distribution  is  actually  from  a 
mixed  and  not  a  true  normal.  It  also  indicates  that  if  the 
variance  is  extremely  large  then  the  affect  of  the  mixing  is 
greater  than  for  the  smaller  values  and  that  only  a  few 
extreme  values  are  required  to  greatly  reduce  the  confidence 
interval  for  the  estimate. 
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APPENDIX  A 

In  order  to  derive  the  equation  for  the  second  order 
approximation  for  variance  of  SEP  (equation  (1))  we  start  by 
noting  that  the  second  order  Taylor  Series  Expression  for  the 
SEP  function  is  given  by  (ignoring  the  remainder  term) : 

SEP(X)  =  SEP(\)+[Diii-LSEP(k)]TAk  +  —  AA.7[r£2S£'^A)]AX 

where  X,  D^^SEP,   and  D^^SEP   are  defined  in  the  text,  X    is  the 

vector  consisting  of  the  maximum  likelihood  estimators  for 
each  element  of  X  and  AX  =  X-X  .   Our  goal  then  is  to  find  an 

approximation  for  the  variance  of  SEP(X).   Because  SEP(X)  is 
asymptotically  unbiased  and  the  mean  square  error 

E[(SEP{X)-SEP(X))2]   =  Var(SEP(X))+(E[SEP(X)]-SEP(X))2 , 
we  can  get  this  approximation  by  looking  at  E[SEP(X)  -SEP(X)  ] 2  . 


26 


Applying  the  tensor  identity  ABC  =  A®Cr£  to  the  last  term 

and  then  transposing  we  can  write  the  second  order  Taylor 
equation  as: 

SEP(X)-SEP(X)   *   [D|t _ j. SEPjX)] rA X  + 1  ( A X ^A X ^ffi, £ SEPjX)] 
=  pB  jjSEF(X)] rA  A,  +  — \DJ_ s SEPjX)] r(A X®  A X) 


Now  we  note  that  since  E(X)  =  X  it  follows  that  £(AX)  =  0  . 
If  we  square  each  side  and  take  the  expected  value  we  see 
that: 


e[sep{X)-sei\X)]2  «  p|ir2.«5£^[A)}3!E[AXAXs]p|ll(2»s,£P[X)] 

qf  T,uggi^A.)]  r£[(A A®  A A.)(A X0A  X)  *][£>,?  T_USEP(X)} 


JLh2 

4 


In  order  to  evaluate  £[(AA.<8>AA.)(AA®AX)7]  we  note  that 

var(AX®AX)  =  2(SR)   (P®P)  (S7?)T  [Ref.  7:p.  44]  where  the  matrix 
SR  has  the  following  qualities  [Ref.  7:  pp.  32-38]: 

SR  is  symmetric 

-  (SR) (SR)  =  SR 

-  2  (SR)(l&tf   =  u®x  +  M®U 
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for  example,  if  u  and  v  are  of  dimension  2  then 


SR  = 


10  0     0 

1  1 

0    —  —    0 

2  2 

1  1 

0    —  —    0 

2  2 

0     0  0     1 


and 


Then   since, 


^(AX^AXXAX^AX)7]  =  var(b\<8>A\)+E{Ak®A\)[E(A\®A\)]T 


AXAXT  =  AX0AX 


if  we   define   P  to  be    i^AXAA.7]  =  var(AA)    we   see   that 

var(AA®AX)+£"(AX®AX)[£(A;i®AA)]T  =  2Sf?(PS>F)(£i?)r+££l 

=  2(SR)T(F®P)(SR)+EEl. 
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Since  D^^SEP     is  a  symmetric  matrix  it  follows  that 
(SR)Du,xSEP  =  DJ^SEP   and  we  see  that: 


E[SEP(X.)-SEP(X))2  *  [D^xSEF(k)]TI[D^zSEP(\)] 

+  \[Dt.  zSEPjk)}  TF®P[Dl  ySEPQ.)] 


which  in  turn  is  approximately  equal  to  the  variance  of 
SEP(X)  . 

Finally,  since  S  is  a  3x3  symmetric  matrix,  we  need  only 
describe  the  unique  terms  (which  we  have  denoted  by  Eu)  to 
arrive  at  the  second  order  equation  for  SEP  given  by  equation 
(1). 
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APPENDIX  B 

A.   APPROACH 

It  is  seen  that  each  of  the  12  Wi    functions  consists  of 

a  double  integral  of  cos  and  sin  functions  multiplied  by  the 
function  g(i\!sin0cos<l>,  -RsinesincJ),  J?cos6) .  Each  of  the  required 
derivatives  carries  through  the  integrals,  so  that  we  need  to 
evaluate  each  in  respect  to  the  function  g. 

We  will  first  evaluate  each  of  the  partials  of  g,  and  then 
carry  these  through  to  the  functions  Wi . 

Recall, 

g(R,Q,$)  =   exp/-— [(i?sin8cos(()-^1)2o11 

+  2(i?sin8cos4)-jjt1)(i?sin0sin<|)-n2)o12 
+  2(i?sin6cos4>-n1)(i?cos8-n3)o13 
+  (i?sin6sin<J)-n2)2o22 
+  2(#sin6sin<t>-n2)(i?cos0-n3)o23 
+  (i?cos6-n3)2o33]}. 
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If  we  define  the  vectors   A, 


A2,   A3    and    A4    by: 


A,   = 


sin28cos24> 
2sin28cos4>sin4> 
2sin0cos8cos(J) 

sin26sin2<|> 

2sin6cos8sin4> 

cos20 


sin0cos<t>m 
sin8sin4>ji1+sin8cos(t)n2 

cosSiij^+sinecos^jia 

sin8sin<t)|ji2 

cos0^2+sin6sin<J)^3 

cos8p,3 


^3  =  [l*i   2^i^2   2^i^3    V-l   2^2H3    M-l]7". 


^4     = 


■RsinScos^-H-L 

i?sin8sin<t>-ji2 

i?cos8-n3 
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Then  treating  R  as  an  independent  parameter  and  taking  the 
partial  with  respect  to  R  we  have: 

a  =  -g'{a11(J?sin26cos2<|)-ji1sin8cos(|)) 

+O12(2i?sin28cos4)sin(f)-jji1sin0sin<t>-jji2sin0cos<|)) 

+a13(2i?sin8cos8cos<|)-p,1cos0-n3sinecos<t)) 

+o22(^sin28sin2<J>-n2sin8sin4>) 

+a23(2i?sin8sin<t)cos8-p.2cos8-^3sin8sin(t)) 

+o33(i?cos28-n3cos8)} 


Similarly,  we  derive: 

-^-   =  gr-{a11(i?sin8cos<j)-ji1)+o12(i?sin8sin<t>-|ji2)+a13(i?cos8-H3)} 
1  =  ^-{[o11   a12   a13]A4} 


3-2-   =  gr-{a12(i?sin8cos<|)-n1)+o22(i?sin8sin(t)-|x2)+a23(i?cos8-H3)} 
'  =  g-{[o12   o22   a23]A4} 


-3-^   =  g-{o13(i?sin8cos4>-^1)+a23(i?sin8sin<t>-n2)+a33(i?cos8-n3)} 
=  gfa13   a23   o33K} 
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In  order  to  derive  the  partials  with  respect  to  sigma 


we 


will  first  generate  the  partials  -3^—  and  the  partials  of — 

dakl  Q 


with  respect  to  akl.      These  values  are  given  below. 


dalx 


do33 


oil 

oil 

oil 

a22 

-a^o11 

-o12axl 

-a^o11 

a33_a220ll 

Q2 

g  -2o1:la12             g33-2a12o12  -^-2a13o12  -2o22o12 

da12  q2  q2 

g  -2o"o13           ^»-2o"o13  ^-2o13o13  ^i3-2a22o: 

ao13                                          g2  g2  g2 

5  a21      iia22             -a12a22  "°13-o13a22  -a22a22 


2 


3a22  g2  g 

2o21-2o11o23     -^i3-2a12o23      -^l2-2o13a23  -2o22a23 


ao23  g2  g2  g2 


a22_all 

a33 

Gl2-a12o 

33 

-a13o33 

an 

-a22o33 

g2 

g2 

sll 

g2 

3 

"<723-a23o11 

g2 

°22_033all 

g2 

-a11 

dolx 

2g 

d 

°13-2o23o12 

g2 

'2°12-2a33a12 

g2 

-a12 

^°12 

g 

a 

°12-2o23a13 

g2 

-2o33o13 

-a13 

5°13 

g 

a 

_a23a22 

°H_a33a22 

g2 

-a22 

5a22 

2g 

-a  -a23 


il-2o23o23  -2o33a23 


do„  a2  g 


'23         g 

_3_  -a23033  -033033  -^ 

aa33  2g 
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We  can  now  use  these  values  to  find  the  desired  partials 


dg 


dalx 


12/rH 


,13_11 


=  g-S—^RsinQcos^-^^ia11) 

+2(i?sin0cos<t>-ki1)(i?sin6sin4>-|i2)o12o 
+2(JRsin8cos4>-n1)(^cos8-n3)a13o 

+(i?sin6sin<t>-M.2)2  a22a 
+2(i?sin8sin<t>-^2)(i?cos8-^3)  a23©11* 


Z21 


23„lla.  "23 


+(i?cos8-n3)2  o33o 


__^22 


rll 


=  ^•p[o11i?2(S-1)uA1-2i?a11(S-1)uA2+o11(S-1)uA3] 


+ 

0 

o  o    "a33 

Q2 

a23 

g2 

_a22 

g2  J 

[R2A1-2RA2+A 

H 

fl 

o11(S"1)u+ 

0    0    0 

"°3 

Q2 

3      °23      -°22Xj 

g2      g2  JJ 

[tf2Ar2#A2+A3]-° 


li 


Following  this   same  methodology   it   is  easy  to   see  that 


dg    _ 

da 


12 


=  g-\o 


12(S_1)U+ 


q         °33         "°23      q        "q!3        °12 

2g2      2g2  2g2      q2 ). 


[i?2A1-2i?A2+A3]-o 


12' 


dg 

ao13 


=  g\ 


o13^"1)"* 


q        ~a23         °22         °13        -°12      0 

2g2      2g2      g2      2g2 


[i?2A1-2i?A2+A3]-o 


13' 


do22        y{2 


a22(S  1)u+|     a33 


0    —A3    o    0    5i 

g2  g2  )\ 


[i?2A1-2i?A2+A3]- 


r22 
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dg    _  „ 

da 


23 


l« 


=  do^S"1) 


°23        -°13        _012      q      _^11      q 


gr2       2q2      2q2  2q: 


[R2Ax-2RA2+A^-a 


23' 


JSL   =  Jl 

3o33        y    2 


a33(S-l)U+|_^22      _^12      0      _^ 

g2       g2  g' 


ii 


0    0 


[iPAi-lBAt+A^S. 
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We  are  now  able  to  calculate  each  of  the  partials  of  theJVj 

by  simply  computing  the  effects  of  the  four  A  vectors  on  the 
original  W  functions  and  applying  the  formulas  derived  on  the 
previous  pages.  We  will  carry  through  the  complete  derivation 
for  Wx   and  then  show  the  solutions  for  the  other  cases.   Since 

the  vector  A3  does  not  involve  trigonometric  functions,  it 

will  not  be  necessary  to  compute  any  partials  for  it. 
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B.   PARTIALS  OF  Wl 

JV,   =  CCX  11!=/      /    cosQcos<bg(R,Qt$)clBd<b    so  we  can  see  that 

'  '  '    Jo    Jo 


the  integrals  of  gA  will  yield: 


I      \     cos0cos4><7*A,  =  /  f  g' 

Jo    Jo  Jo    Jo 


sin20cos8cos3(|) 
2sin28cos0cos2<J)sin4) 

2sin8cos2Gcos2<j> 

sin26cos0cos(|)sin2<t) 

2sin8cos26sin<t>cos(j) 

cos38cos<|) 


cc  -cc 

'■''■'1,1,3,1     ^^3,1,3,1 


*(CSlt  1,1,1     CS3f  1,1,1     CSi,i,3,i+CS3,i#3#l) 

^(•^<-^l,l,2/l~'^^'3,l/2,l) 

CC  -CC  -CC  +CC 

'-'■'1,1,1,1     '-'-3, 1,1,1     *-*-!, 1,3,1     ^^3,1,3,1 


•^1,1,1,2     ^3,1,1,2 


CC 


3,1,1,1 


/■2it  /•  n  r2it  rn 

/    cos8cos(j>g-A>  =1  g 

0     J  o  Jo     JO 


sinScosScos2^! 

sin8cos8cos4>sin<|)^1+sin0cos8cos2<t)|i2 

cos28cos(|)p.1+sin8cos8cos2<l)p.3 

sin8cos8cos4>sin<J>n2 

cos28cos(J)|i2+sin8cos8cos(t)sin(l)p,3 

cos20cos<t>M.3 

^l^l,2,l,2+2h^l,2,2,l 

AM^.i.i.i+^s^i.a^.i 

1*2^1,2,1,2 

4^2^2,1,1,1+^3^1,2,1,2 
4^3^2,1,1,1 
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I      /    cosecos^g'A.   =  /      /    <7* 
Jo    Jo  Jo    Jo 


i?sin8cos6cos2<|)-cos8cos(|)^1 
.Rsin8cos8cos<t>sin4>-cos8cos<}>p.2 

i?COS20COS<|>-COS0COS<|)H3 


—  sin28cos2(|)-ji;Lcos8cos(t) 

—  sin28sin2(()-|i2cos8cos(J) 
-Rcos28cos<t>-ii3cos8cos4> 


^l,2.1.2^-4^2^1,l,l.l 


For  future  use,  we  point  out  that  the  partials  of  A2   andA4 

involve  only  3  distinct  trigonometric  terms,  which  we  will 
define  by: 


B1  =   sin6cos4> 
B   =  <JB2   =  sin6sin<|> 
S,  =  cos8 


then  we  see  that 


A,   = 


B2\i1+B1\L2 
B3[i1+B1\k3 

B2»2 

B3\i2+B2\i3 

B^2 
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and 


^4     = 


BxR-\ix 
B2R~\i2 
B3i?-n3 


Thus,  we  need  only  compute  the  effects  of  Ax  and  B  on  each 
of  the  12  W  functions  to  obtain  the  desired  partials. 

C.   PARTIALS  OF  W2 


W2   = 


/    cos2>Qcos<bg(R,Q,4>)dQd<b   so  the  partials  are 

o    Jo 


given  by: 


cossecos^^  = 


sin28cos38cos3<t>    . 

2sin28cos38cos2<|>sin(J) 

2sin6cos8cos36cos2(t> 

sin26cos36sin24>cos(J> 

2sin8cos6cos38sin<|)cos<|) 

cos28cos38cos<J) 


2  ^-"1 ,3,3,1" C*-*l ,  5 ,  3 , 1  ~  * C i ,1,3,1 
2 (2  CSX  #  3  #  1#  j_ ~  CS^  ,  5/  j.,  j.  ~  C^  /  ! ,  ! ,  !  ~2  CS^  ,  3 ,  3 ,  !  +  CS^  #  5 ,  3  #  !  +  CS1 ,1,3,1) 

2('5<-l,  5,  2, 1_,5C1, 1,  2,  l) 
^  ^-1 ,  3  , 1 , 1  ~  *~  C 1 ,  5 , 1 , 1  ~  ^-  *- 1 , 1 , 1 , 1  ~  ^  ^~1 ,  3 ,  3 , 1  +  ^l ,  5 ,  3 , 1  +  ^-1 ,1,3,1 

■^1, 5,1, 2  ~&Sl,  1,1,  2 

rr  +2pr  +rr 

'-'-1,5,1,1  T «■«-«-!,  3,1,1     t-t-l,  1,1,1 
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cos38cos(J)B  = 


sin6cos38cos24> 

sin6cos38cos4>sin<t> 

cos6cos38cos<t> 


2(-SHl,4,2,1     ^-l,2,2,l) 

S3  -SS 

2(cc1/4#1#1+cq/2#1/1) 


PARTIALS  OF  W3 


W, 


/'  2  it  /*  w 
I    cos8sin<J)5f(i?/8/<|))d8d(J)    so  the  partials  are 
0    Jo 


given  by: 


cosSsirKf)^  = 


sin28cos8cos2(J)sin(t) 

2sin26cos6cos<J)sin2<J) 

2sin8cos26cos(|)sin<t> 

sin28cos8sin34> 

2sin8cos28sin2<t> 

cos38sin(J) 


CSi ,  1 , 1 , 1  ~  Ci?  3  , 1 , 1 , 1  ~  CSj, ,  1 ,  3  , 1  +  ^ ^3  ,1,3,1 

^  ( ^H. ,  1 , 1 , 1 "  CCX  f  1  #  3  f  1  ~  CC3  ( 1  ( 1  _  1  +  CC3  ,1,3,3.) 
^1,1,  1,2  -"S^a  ,1,1, 2 

^1 ,1,3, 1  ~C^3  ,1,3,1 

2(5'«S1 ,1,2,1  "^3  ,1,2,1) 


OS 


3,1,1,1 


cos8sin(J)B  = 


sin6cos8cos(t>sin<|) 

sin8cos8sin2<J> 

cos28sin(j) 


SS. 


1,2,1,2 


2SS. 


1,2,2,1 


4  OS 


2,1,1,1 
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E.   PARTIALS  OF  W4 


WA 


/    cos3dsin<bg(R,  6,  $)dQd<\>      so    the    partials 

o    Jo 


are   given  by: 


cosBOsin^AL  = 


sin20cos30cos2(t>sin(|> 

2sin20cos36cos(|>sin2(t> 

2sin8cos0cos30cos(|)sin<J) 

sin20cos30sin3<|> 

2sin0cos0cos30sin2<J) 

cos28cos30sin<t> 


^  *- "^l ,  3 , 1 , 1  ~ 2  ^  1 ,  3 ,  3 , 1  ~  ^  1 ,  5 , 1 , 1  +  *- ^  1 ,5,3,1"  C ^1 , 1 , 1 , 1  +  "*1 ,1,3,1 

2  (2  CC^ ,  3 , !  7 !  ~  2  C^  /  3  ( 3  ( !  ~  CCX  >5>ltl  +  CC1  r5,3,i~  CC1  tliltl  +  CC1  i1i3i1) 

2  C°\ ,3,3,1"  ^  °1 ,  5 ,  3 , 1  ~  ^1 ,1,3,1 
^  ( ^1 ,  5 ,  2 , 1  ~  ^1 , 1 ,  2 , 1 ) 

^1,5,1,1+2^1,3,1,1  +  ^1,1,1,1 


cos30sin<t>B  = 


sin0cos30cos(t>sin4> 
sin0cos30sin24> 
cos0cos30sin<j> 


S&1,4, 1,2  ~^1,2,1,2 
2(5^ ,  4  ,  2 , 1  ■"•S'S'i ,  2 ,  2 ,  l) 

2(C5X  4 #1  ^i+CS^ /2/i,i) 
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F.   PARTIALS  OF  W5 


W5   = 


I  sinQg{R,  0,<|>)d0d(t>  so   the  partials   are 

o    jo 


given  by: 


sinS^  = 


sin30cos2<J> 
2sin30cos<t)sin<t) 
2sin26cos0cos<J) 

sin30sin24> 

2sin26cos8sin<|) 

cos20sin8 


*^3  ,1.2,1 
S3  3  ,1,1,2 

2(CC1  r  x,  i  ,  i  ~  CC3  r  x ,  i ,  i) 

°&3  ,1,2,1 

2(C"Sx  r  x,  X,  1  ~ ^3  , 1 , 1 ,  l) 

'-'^1 , 1 ,  0 ,  0  ~^-3  , 1 ,  0  ,  0 


sin0B  = 


sin20cos(}> 
sin20sin<|) 
sin0cos0 


25C 


2,1,1,1 


255. 


2,1,1,1 


sc. 


1,2,0,0 
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PARTIALS  OF  W6 


^6     " 


/    sin2Qg(R,6,$)dQd<\>    so     the     partials     are 

o    Jo 


given  by: 


sin20A1  = 


sin26sin20cos2<t> 
2sin28sin28cos(J)sin(t) 
2sin0sin26cos6cos(l) 

sin26sin26sin2<J) 

2sin0sin29cos0sin<|> 

cos20sin20 


^  ^-1 ,2,2,1  ~SCl  ,4,2,1 

2SSlt  2, 1,  2~^X,  4, 1,  2 

2(^1,  2,  1,1~^1,  4, 1,1 +2  ^^2,1,1,1) 

2  «!  #  2  /  2 , 1  ~  '5'S'i ,  4 ,  2 , 1 

2(C^l,  2, 1, 1  ~CSl,  4, 1, 1  +  25.S2> 1, 1,  l) 

^  ^-1 ,2,0,0  +^-i  ,4,0,0 


sin20S  = 


sin0sin20cos4> 

sin0sin20siri(t) 

sin20cos0 


CC  -CC 

ot-i,  1,1,1    ^-^-1,3,1,1 

t-,Ji,  1,1,1    *-*Jl,3,l,l 

'-"■'1,3,0,0     '"'^l, 1,0,0 
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H. 


PARTIALS  OF  W7 


I    s in3 Qg(R,  6  ,ty)dQd<b    so     the     partials     aire 

o    Jo 


given  by: 


sin36A1  = 


sin20sin30cos2<f> 
2sin20sin38cos<l)sin<t> 
2sin8sin30cos9cos(J) 

sin20sin30sin2<|> 

2sin0sin30cos0sin<t> 

cos20sin30 


2  S^l  ,3,2,1  ~SC-i  ,1,2,1  ~SCl  ,5,2,1 

^  (  CC-L ,  1 , 1 , 1  ~  ^"1 ,  5 ,  1 , 1 ) 
^  ^^1 ,  3  ,  2 , 1  ~  ^1 , 1 ,  2 , 1  ~  ^1 ,  5 ,  2 , 1 

2(^1,  1, 1,  \~^^X,  5, 1,  l) 
o(-l,5,0,0  ^oul,  3,0,0  OL-l,  1,0,0 


sin30S  = 


sin0sin30cos(f> 

sin0sin30sin4> 

sin30cos0 


ut-l,2,l,l  0<-l,4,l,l 
^ "^1 ,  2 , 1 , 1  ~  ^ ^  1 ,4,1,1 
'^^-l,4,0,0+'^<--l,2,0,0 
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PARTIALS  OF  W8 


W0 


/*  2tc  /*  it 
j    sin2Qcos<bg(R,Q,$)dBd<b   so  the  partials  are 
o    Jo 


given  by: 


sin^cos^AL  = 


sin48cos34) 

2sin48cos24>sin<t> 

2sin38cos6cos2<t> 

sin40sin24>cos<J) 

2sin36cos8cos<|>sin<t> 

cos28sin28cos<t> 


QSC4, 1,3,1 

16(SS4#1  /1#1-554  /1/3/1) 
2(2 5CX  2  2 tl~SCx  4  2  x) 

®  (SC* ,  1 , 1 , 1  ~  °^A  , 1 ,  3 ,  l) 

"n,2,l,2~^l,4,l(2 
^(•^^2 , 1 , 1 , 1  '3^4  , 1 , 1 ,  l) 


sin26cos<f>.B  = 


sin38cos24> 
sin36sin<t>cos4> 
sin28cos8cos<t> 


2SC 


3,1,2,1 
^3,1,1,2 
2(CC1  f  1;  1#  x~  C^-3  , 1 , 1 ,  l) 
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J.       PARTI ALS    OF   W9 


Wa 


/    sin2Qcos2^>g{R,Qf^)dQd<\>   so  the  partials  are 

o    jo 


given  by: 


sirr'Ocos2^!  = 


sin50cos44> 
2sin50cos34)sin<t> 
2sin40cos0cos3<t> 
sin50sin2<|)cos2<t) 
2sin40cos0cos2(J)sin<t) 
cos20sin30cos2<J> 


4^-5,1,4,1 

SSSi  1#  1#  4  +2SS5i  lt  Xi  2 

°(<-^l ,  1, 3, 1  "^  CC3  (  j  #  3  <  j  +  CC5 f  i,  3,  i) 

4(*^5 ,1,2,1  ""S^s , 1 ,  4  ,  l) 

°{CS\ ,  1 , 1 , 1  ~  ^Sj, ,  1 ,  3  , 1  ~2  Co3 , 1 , 1 , 1  +2  CS3 ,1,3,1  +  CS^  f  1#  1#  j  ~  CS5  ,1,3,1) 

*(^^3,l,2,l-^H5,l,2,l) 


sin30cos2<J)S  = 


sin40cos34> 
sin40sin4>cos2<J) 
sin30cos0cos2(t>. 


85C4, 1,3,1 

8  {SS4 , 1 , 1 , 1  ~  ^4 ,1,3,1) 

<il-'v-'l,2,2,l     "-"-1,4,2,1 
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K.       PARTI ALS    OF   W10 


W. 


10 


=  SS2  1  1  1  -  I      \    sin20sin(j>sr(i?,0 ,  <f>)d0d(J>   so  the  partials  are 

Jo    Jo 


given   by: 


sirr^SsiiHf)^  = 


sin40cos2<j>sin<t> 

2sin48cos4>sin2(J) 

2sin30cos0cos<t>sin.(J> 

sin48sin34> 

2sin30cos0sin24> 

cos20sin20sin<}) 


8(554  ,1,1,1  ~SSa  ,1,3,1) 
16(SC41  #1  tlSC4  13  tl) 

2  S&i ,  2 , 1 , 2  ~  ^1 ,  4 , 1 , 2 


855 


4,1,3,1 


^  S^i  ,2,2,1     ^  55x  ,4/3,1 
8(^2. 1, 1,1  ~SSi,l,  1,  l) 


sin20sin4>B  = 


sin30cos<t>sin(t> 

sin30sin24> 
sin20cos0sin<J> 


^3,1,1,2 

1 

2 

^  "3  ,1,2,1 

2(C51(1  1(1-C531#1  x) 
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L.   PARTI ALS  OF  Wll 


J    sin36sin2(t>g<i?,e/<|))d0d<t)    so    the    partials 

o    Jo 


are   given  by: 


sin30sin2<|>.A1  = 


sin50cos2<t>sin2<t> 
2sin50cos<t)sin4>sin2<t) 
2sin40cos0cos(J)sin2<J) 

sin50sin24>sin2(J) 
2sin40cos0sin<t>sin2(j> 

cos20sin30sin2<|) 


•^5 , 1 , 1 ,  4  +  ^  •£•£ 5 ,1,1,2 


2(2SC512tl    5C51jl4    ^--5,1,1,2) 

^l,l,l,3  +  ^l,l,l,r^('^3,l,l,3"^^3,l,l,l+^5,l,l,3+('^5,l,l,l) 


2^5,1,1, 2     ^5,1,1,4 


4(^-1, 1,1,1     ^-"1,1,1,3     ^^-<-3,l,l/l+^^'-3/l/l/3  +  ^-^5/l,l/l     ^-^-5,1,1, 3) 

4(^3,  1, 1,  2  ~SSs,  1, 1,2) 


sin30sin2<J>jB  = 


sin40cos(t>sin2(|> 
sin40sin4>sin2({> 
sin30cos0sin2<}) 


4(^4. 1,1. 3+^4, 1.1,1) 

^(■^Q ,  1 , 1 . 1  "^^ ,  1 , 1 , 3) 
2  S*i ,  2 , 1 , 2  "  ^1 , 4 , 1 , 2 
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M.   PARTI ALS  OF  W12 


/    sin30sin24>3'(.R/0/<|>)d0d<t>    so    the    partials 

o    Jo 


are   given  by: 


sirr'Ssin^A,^  = 


sin56cos2<l>sin2<j) 
2sin58cos<t>sin3<i> 
2sin40cos0cos<|>sin2(t> 
sin50sin44> 
2sin40cos0sin3<t> 
cos28sin30sin24> 


^(^®5,l,2#l-"^5#  l,4,l) 
2 SSg,  1,1,  2  "^5,1,1,4 

°\pC\  ,1,1,1""  CC-i ,  1, 3, 1  ~2  CC3  ilflil+2  CCj  f  1 1 3  ( 1  +  CC5  tlilil~  CC5  ilt3il) 

^  ^5, 1,4,1 

°  (  ^1  ,  1 , 3 , 1  ~  ^  CS3  /  x  ^  3  f  x  +  CS5  i1i3i1) 

M        3. 1.  2,  1_,^5,  1,  2,  l) 


sin30sin2<J)B  = 


sin40cos4>sin2(f> 

sin40sin3<J) 
sin30cos0sin24> 


°(^-4,l,l,l-^-4, 1,3,1) 

QSS4, 1,3,1 
^  -SOl ,  2  ,  2  , 1 "  ^1 ,4,2,1 
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APPENDIX  C 

A.   PARTIALS  OF  U1,U2  AND  U3 

This  step  simply  involves  taking  derivatives  of 


T    - 


[Ul     U2     U3J        = 


W8        V10  W6 


W^ 


IT 


w5      2w5 


with  respect  to  each  parameter.   The  results  are: 


dux 
dw~ 


w, 


w. 


du 


1  _   1 


dwe         w5 


du 


2      _  V10 


dw„ 


«? 


du 


2       _ 


^10 


w* 


du 


3     _ 


Wc 


dwc 


2w: 


du 


3      _ 


dw,         2  wc 


and  all   other  partials  =   0. 
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B.       PARTIALS    OF   U4,..,U9    WRT   MU 


o11 

3u4 

=    .       "8 

o12 
o13 

d» 

2w5 

o11 

du5 

-   ..  Vi° 

o12 
o13. 

d\x 

W5 

o11 

du6 

-   .     we 

o12 
[a12 

dp 

2ws 

o12 

du7 

-   ..  vio 

o22 
o23 

d\i 

2ws 

o12 

due 

-   .     We 

022 

o23 

dp 

2w5 

o13 

du9 

-    .       V6 

O23 
O33 

d\i 

4ws 
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C.   PARTIALS  OF  U4,..,U9  WRT  SIGMA 

To  solve  for  these  partials  we  use  the  table  of  partials 
given  in  the  previous  section  and  use  standard  calculus  rules: 
1.   Partials  of  U4 


du 


4       _ 


da 


-o11^ 


11 


du 


4       _ 


da 


12 


'33 


fRwxl     p2We\    _    (  R(w1-w2)_  \i,w8) 

a,.,  oT;  °23  Qr.,  ->,., 


4wc       2  wK 


8wc 


2w, 


5    / 


au 


5oi: 


4       _     _ 


<72 


/jg(w1-w2)_  H3*0_a    f  Rvr1±  _  ]i2we' 


22y      8wc  2vc 


4  Wc       2w, 


5  y 


au 


4       _ 


a<? 


22 


a22u4-^ 
Q2 


'13 


/^v^-^)_M^g\_     f  fly9  _  Hj^, 


8wc 


2w, 


5  / 


2\2w5      2w5 


du 


4       _ 


da 


=  -2o"uA- 


23 


4  2 


f  Aur,-^*^ 


'23 


W* 


-a 


13 


f  g^_    H2*M_„      f  jjfrV^g).   ^M 


4  w5       2w5 


-a 


12 


8wc 


2w, 


s  / 


au 


4       _     _„33 


da 


=  -oJJu,-. 


33 


'12 


f^ll_hV8^_0      ^^9_^1V8 


4w5       2w5 


'22 


2w5      2w5 


51 


2.   Partials  of  U5 


du 


5       _ 


a°ll 


-O11^ 


du 


5      _    _o«12 


da 


=  -2o12uc- 


12 


5  2 


'33 


°V12  _   ^2V10  | 

°23 


V      ^5 


W= 


^3-^4)  _  m^ 


4w~  w, 


s     /J 


du5 
6o13 

-2o13u5-^- 
<72 

du5 

3°22 

-022u       _L 

<72 

'22 


^^3-^4)_  ^3^10^_a      (  RW12_  \l2W1Q) 


4w*  Wk 


'23 


wc  w, 


5    ; 


'13 


Ar.,  ,.,  033  Or,  TT 


4WC  W, 


5    y 


y     2W5  W5        j 


du 


5       _ 


5a23 


-2o23u5- i- 
Q1 


'23 


^^n  "2  l^i^io^ 


wc 


-a 


13 


fjgg^2-jW0>L,T      f  jgVj^),    ^3^10 


VVC 


-a 


12 


4wc  wc 


du 


L     =     -«33JT    _ 


3o 


=  -oJJu 


33 


5  2 

<?2 


'12 


*V12         ^2^10^      „       f^Vll         HlVl 


Wc  V. 


-o 


22 


5    y 


2^        w^ 


3.      Partials   of  U6 


du 


6        _ 


3o 


-o11^ 


11 


du 


«       =     -Orrl2 


3o 


=  -2o12u£- 


12 


6        2 
<72 


'33 


fjggV2^1_    1*2*6  \_a      I  RiW7+W5)_   ^3^6 


4  v.-  2  Wc 


23 


4Vc  2  Wc 
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du 


6       _     _->„13 


da 


=  -2o1Ju£- 


13 


6  2 


(  R(w1+w5)_  y^wA         f  R(w3-w4)_  \i2we) 

°22  A...  ->,.,  °23  7TT  o.. 


4vq  2wc 


4  Wc  2  W, 


5  yj 


au 


6       _ 


da 


•o22u6— i- 


22 


"13         a...  o,.,         °33         777  ~oT7" 


4(/c  2wc 


4  v/c  2vc 


du 


6       _ 


3o 


23 


-2a^U6-^L 


'23 


^l-^2)  _   ^1^6 


2wc 


wK 


-a 


13 


f  Jfofr-Wfr).  1*2*6^ _,      (R{W1+W5)_   ^3^6 


4wc  2vc 


12 


4wc  2wc 


au 


6 


ao 


33 


'12 


5 


-4)_  ^ 


w,-w. 


2*0_„      f  j&V^),   ^1*6 


4W.  2  Wc 


-o 


22 


4  wc  2  wc 


4.   Partials  of  U7 


du7 


da±1 


-«"m    -. 


=  -o^^u. 


<T 


'23 


f^3-^4)^3^10^_    (Rw12_  \i2wl0 


8vc  2  w, 


'33 


s   / 


v  2w5        2w5 


du7 


=  -2o"u 


,  -    1 


<r 


'33 


f  Rw1±_  ^i^iq>|_ct    p(^3-^4)_  H3*V 

/I    t.t  O  r.r  13  o  T.T  O  T.» 


V  4^ 


2w, 


5     / 


8w,  2  w. 


5     / 


au7 


5oi: 


-2o"u7-^- 


/i?v^12-ji2^10\ 


'13 


Wc 


-a 


23 


■^11      M'l^io 


8  w7  2w5 


53 


du7 


■o22Hj 


du7 


do 


23 


-2a23u7— ±- 
q2 


'11 


R(w3-Wi)    \x3w10\        (R*ii_  ^1^10 ' 

O  r.r  1  T.r  13  A   ,.,  O  ... 


8  WK  2wc 


4wc        2  w, 


V    "5 


5  yj 


du7 


3o 


33 


RW11_  1*1*0  [jg^2_jWL0 

4wc        2w,  X1     2w,        2v, 


5.      Partials   of  U8 


d"8 
doxl 


Q2 


'23 


R[w^w5)  _  1*3*6  ^_„    f  j(jV^)_  k*2v^6 


4w5  2  v, 


-o 


33 


5    ) 


4  VVB  2v, 


3u 
3^ 


8       =     _0«12 


=  -2o12uc 


12 


'33 


fjggCggl,   ^1^6  \_Q      (  R(W1+W5)_   ^3^6 


4  wc  2Wc 


13 


4  wc  2  w, 


_3u 


8       _     _0/rl3 


=  -2ouu0- 


13 


'13 


R(w3-wA)_  \x2w6 

2Wr  We. 


(R{w1-w2)_  \^wA        (Rjw7+w5)_  \i2w6 

°23  TT,  ~K7T\     °12  TTTr  <->,., 


4wc 


2ws  )       12[       4W5 


2K 


du 


8       _ 


do 


-o22uc 


22 


8<j 


8       _ 


■2o23uQ-. 


23 


11 


RQ*,+w5)_We)         (Riwi-W2)_  ^i^6 

/I    T.r  O  r.r  13  /I   ...  O  ... 


4  v-  2  w, 


5   / 


4  v/c  2  vc 
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6Ur 


do 


■o33u0-A 


33 


'12 


R(w1-w2)_  1*1*0         fj^V3)_jV^x 

A   r.r  1  r.r  11  A   t.r  O  r.r 


4wc  2  w, 


5    / 


4wK  2w, 


5    /J 


6.   Partials  of  U9 


du^ 


da 


■o11^-^- 


ii 


'23 


ig(V3-V4)_   We)-      (R(W7+W5)_   1*3^6 


8wR  4  wc 


'22 


8  V-  4  WK 


du9 

3°12 

<?2 

°12 

4V5 

^3^6 

2w5) 

"°13 

r^V3-^4)_ 

{       8ws 

4w5  j 

"°23 

1*1*0" 

4V5   )\ 

du9 
3o13 

-2o-u9-^ 
<72 

°22 

(R(w1-w2)_ 

{          8W5 

1*1  V 

4V5    , 

|-°12 

I         8V5 

\*2w6  y 

4V5   JJ 

du9 

3°22 

aJ^l-N- 

"I       8*5 

4^5  J 

8^5 

**3»0 

^5    jj 

du9 

5°23 

-2o23u9-A 

aii 

(R(W3-W4) 

{      8"s 

1*2^6 

4W5 

\ 

-°12 

/ 

fB(wl-wi)_ 

{          8W5 

l*i  *0" 

4w5)\ 

dus 


da 


=  -o33u, 


33 
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D.       PARTIALS    OF   U4,..,U9   WRT   R 

du4  i 


dR         8w5 


(4o11w9+2o12w11+o13(w1-w2)) 


-s*  =  TTr(2allvn+4al2vi2+°13(v3-v4)) 


dR  4w5 


1JR      =    ^(al\Wl-W2)  +  al\W>-W*)+0l\W->+W5)) 


2 

"3r  =  -i^(al2^"V2)+°22^3"V4)+°23^7+V5)) 


2. 


E.   PARTIALS  OF  U4f..,U9  WRT  W 
1.   Partials  of  U4 


du4  _  Ro 


13 


dw1  8  w5 


du4    =  _  i?o13 
dw2  8  w5 


56 


a*5       ^5  4 


£  ■  -2fe(oii^+oi2^+o13^ 


'5 


du4    =    flo11 
dw9  2  w5 


d"4     _    i?a12 


dw^  Aw5 


2.   Partials  of  U5 


du5    =    £o13 

dw3  4  w5 


du5     =    _Ra13 

dw4  4  ws 

*±  =  -Au 

dw5  w5    5 

dw13L         2  w5 


57 


d"5     m    Ra12 
dw12  w5 


3.      Partials   of  U6 


du6    _    Raii 


dw±  4  w5 


du£  Pa11 


6    _      Ra: 


dw2  4  w5 


du6    =    Rq12 
dw3  4  w5 


dllc  Pn12 


6    _      Ra- 


dw4  4  w5 


duc  1  p„13 


£   =  -J.Mt»L 


dw5  v5    6      4w5 


||   -  -^(."Hx—I^m) 


du6    .    Ra12 
dw7  4w5 
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4.   Partials  of  U7 


Buj    =    Ro22 
dw3  8  w5 


3^7  Ro 


23 


dw4  8  w5 


a*5        ^5  7 


■fe  ■  --^k^-22^2^) 


du7     _    2?o 


12 


3vxl  4w5 


5u7  i?a 


22 


a^12      2w5 


5.   Partials  of  U8 


8wx         4  w5 


dw2  4  w5 


dw3  4  w5 
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d"e    =  _Ro22 
dw4  4  w5 


duQ  1  Pa23 


£  =  -AUo  +  **: 


d^5  V5      8        4V5 


du 


S  =  -^-(a«n1+a"|ia+o»|i3) 


dw6  2w5 


Partials  of  U9 


du8    =    Ra22 
dw1  4  w5 


du9 

Ra13 

dwt 

8w5 

du9 

Ra12 

dw2 

8w5 

du9  Ro22 


dw3  8  w5 


du9     =    _  J?q23 

dw4  8  w5 


aU9       -    _    1    „    +  *°33 


3^5  V5      9         8u% 


60 


duc 


=  -TTr(°13^+(j23^+o33^) 


dwe  4  w5 


du9  2?033 


dw7  8  w5 
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APPENDIX  D 


PROGRAM  SEPCMP 

PROGRAMMER:    LT  ARTHUR  F.  BROCK 

DATE:  12  MAY  1991 

LAST  MODIFIED:  29  AUGUST  1991 
PROGRAM  DETERMINES  FIRST  AND  SECOND  ORDER  APPROXIMATIONS  TO 
VARIANCE  OF  SEP  GIVEN  VALUES  OF  MU  AND  SIGMA  (VARIANCE-COVARIANCE 
MATRIX)  AS  INPUTS.   OTHER  INPUTS  ARE  PARAMETERS  FOR  USE  IN  THE 
IMSL  SUBROUTINE  QTWODQ  USED  TO  PERFORM  GAUSSIAN  QUADRATURE 
EVALUATION  OF  DOUBLE  INTEGRALS. 

INCLUDE  'COM  DEF' 
INCLUDE  'WVEC  DEF1 
INCLUDE  'PMUCOM  DEF' 
INTEGER  I,OT1,OT2,SAMPSZ,K 

REAL  DT1(10,9),DT2(22,10),DT3(9,22),SEPM(9,9),PMU(3,3),PSIG(6,6) 
C       ,SEPMU(3),SEPSIG(6),CI(2),MPT(3,3) 
CALL  INIT(SAMPSZ.MPT) 
WRITE(16,51)  R,NUMTR 

URITE(16,52)(MU(I),I=1,3) 
DO  20  1=1,3 

URITE(16,53)(MSIG(I,J),J=1,3) 
20  CONTINUE 
OT1=0 
OT2=0 

DO  10  I=1,SAMPSZ 
CALL  INIT2(MPT) 
CALL  WVALS 

CALL  DT1VAL(DT1,SEPMU,SEPSIG) 
CALL  DT2VAL(DT2) 
CALL  DT3VAL(DT3) 
CALL  SEPEST(DT1,DT2,DT3,SEPM) 
CALL  PMUS(MSIG,PMU,PSIG) 

CALL  BNDEST(SEPM,PMU,PSIG,SEPMU,SEPSIG,CI) 
IF(SEPSET  .LT.  (R-CI(1))  .OR.  SEPSET  .GT.  (R+CI(1)))  THEN 
OT1=OT1+1 
IF(SEPSET  .LT.(R-CI(2)).OR.  SEPSET  .GT.(R+CI (2)))  THEN 

OT2=OT2+1 
END  IF 
END  IF 

PRINT*, 'I  =  ',I,REAL(I-OT1)/REAL(I),REAL(I-OT2)/REAL(I) 
WRITE(16,50)  R-CI(1),R+CI(1),OT1,R-CI(2),R+CI(2),OT2 
10  CONTINUE 

WRITE(16,55)  SAMPSZ,REAL(OT1)/REAL(SAMPSZ),REAL(OT2)/REAL(SAMPSZ) 

50  F0RMAT(1X,2('   ( ' , F7.2, ' , ' , F7.2, ' ) ' ,2X, 13, 2X)) 

51  FORMAT(1X,'TEST  FOR  SEP  =  ', F7. 2,/, 1X, 'SAMPLES  PER  TRIAL  =  ',12) 

52  F0RMAT(1X,'         MU  ' ,/,3(F7. 2, ','),/,/, 1X, 
C         '  SIGMA') 

53  FORMAT(1X,3(F9.3,3X)) 

55   FORMAT (1X,/,1 SAMPLE  SIZE:  ■ , I4,3X,2(F5.3,8X)) 
STOP 
END 
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SUBROUTINE  INIT(SAMPSZ.MPT) 

INITIALIZE  FILES,  READ  INPUT  DATA  AND  INITIALIZE  VARIABLES. 

INCLUDE  'COM  DEF' 
INCLUDE  'PMUCOM  DEF' 
INTEGER  N,SAMPSZ,RK 
REAL  MPT(3,3),TOL,NUMSIG(3,3) 
EXTERNAL  RNSET,DCHFAC 
OPEN(15,FILE='/MUSIG  DATA  AT) 
OPEN(16,FILE='/OUTDT  DATA  A1 ' ) 
READ(15,*)  MU(1) 
READ(15,*)  MU(2) 
READ(15,*)  MU(3) 
READ(15,*>  SIG(1) 
READ(15,*)  SIGC2) 
READ(15,*)  SIG(3) 
READ(15,*)  SIG(4) 
READ(15,*)  SIG(5) 
READ(15,*)  SIG(6) 
RE AD (15,*)  ERRABS 
READ(15,*)  ERRREL 
READ(15,*)  IRULE 
READ (15,*)  NUMTR 
READ (15,*)  SAMPSZ 
READ(15,*)  ISEED 
READ(15,*)  NUMSIM 

CALL  RNSET(ISEED+INT(MU(1)+SIG(1))) 

DET=SIG(1)*(SIG(4)*SIG(6)-(SIG(5)**2))-SIG(2)*(SIG(2)*SIG(6) 
C       -SIG(3)*SIG(5))+SIG(3)*(SIG(2)*SIG(5)-SIG(3)*SIG(4)) 
IF  (DET  .LE.  0.0000000001)  THEN 
WRITE(16,50)  DET 
STOP 
END  IF 

SIGINV(1)=(SIG(4)*SIG(6)-(SIG(5)**2))/DET 
SIGINV(2)=(SIG(3)*SIG(5)-SIG(2)*SIG(6))/DET 
SIGINV(3)=(SIG(2)*SIG(5)-S1G(3)*SIG(4))/DET 
SIGINV(4)=(SIG(1)*SIG(6)-(SIG(3)**2))/DET 
SIGINV(5)=(SIG(3)*SIG(2)-SIG(1)*SIG(5))/DET 
SIGINV(6)=(SIG(1)*SIG(4)-(SIG(2)**2))/DET 
DENOM=(SQRT(2.0*PI)**3)*SQRT(DET) 
MSIG(1,1)=SIG(1) 
MSIG(1,2)=SIG(2) 
MSIG(1,3)=SIG(3) 
MSIG(2,1)=SIG(2) 
MSIG(2,2)=SIG(4) 
MSIG(2,3)=SIG(5) 
MSIG(3,1)=SIG(3) 
MSIG(3,2)=SIG(5) 
MSIG(3,3)=SIG(6) 
CALL  FINDSEP 
MSIG(1,1)=SIG(1)/NUMSIM 
MSIG(1,2)=SIG(2)/NUMSIM 
MSIG(1,3)=SIG(3)/NUMSIM 
MSIG(2,1)=SIG(2)/NUMSIM 
MSIG(2,2)=SIG(4)/NUMSIM 
MSIG(2,3)=SIG(5)/NUMSIM 
MSIG(3,1)=SIG(3)/NUMSIM 
MSIG(3,2)=SIG(5)/NUMSIM 
MSIG(3,3)=SIG(6)/NUMSIM 
TOL=100.0*DMACH(4) 
CALL  DCHFAC(3,MSIG,3,TOL,RK,MPT,3) 
CALL  SETDAT 
CALL  MAKES 
50   FORMATdX,' ILLEGAL  INPUT  MATRIX,  DETERMINANT  =  \F8.5,'  PROGRAM1 
C      '  TERMINATED.1) 
RETURN 
END 
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SUBROUTINE  INIT2(MPT) 

*       INITIALIZE  VALUES  USED  EACH  TIME  THROUGH  LOOP. 

REAL  MPT(3,3),TMU(3),TSIG(6) 
INCLUDE  'COM  DEF' 
INCLUDE  'PMUCOM  DEF' 
20   CONTINUE 

SIG(1)=0.0 

SIG(2)=0.0 

SIG(3)=0.0 

SIG(4)=0.0 

SIG(5)=0.0 

SIG(6)=0.0 

MU(1)=0.0 

MU(2)=0.0 

MU(3)=0.0 

DO  15  I=1,INT(NUMSIM) 

CALL  TRIALS(NUMTR,MPT,TMU,TSIG) 
SIG(1)=SIG(1)+TSIG(1) 
SIG(2)=SIG(2)+TSIG(2) 
SIG(3)=SIG(3)+TSIG(3) 
SIG(4)=SIG(4)+TSIG(4) 
SIG(5)=SIG(5)+TSIG(5) 
SIG(6)=SIG(6)+TSIG(6) 
MU(1)=MU(1)+TMU(1) 
MU(2)=MU(2)+TMU(2) 
MU(3)=MU(3)+TMU(3) 
15      CONTINUE 

MSIG(1,1)=SIG(1) 
MSIG(1,2)=SIG(2) 
MSIG(1,3)=SIG(3) 
MSIG(2,1)=SIG(2) 
MSIG(2,2)=SIG(4) 
MSIG(2,3)=S1G(5) 
MSIG(3,1)=SIG(3) 
MSIG(3,2)=SIG(5) 
MSIG(3,3)=SIG(6) 

DET=SIG(1)*(SIG(4)*SIG(6)-(SIG(5)**2))-SIG(2)*(SIG(2)*SIG(6) 
C         -SIG(3)*SIG(5))+SIG(3)*(SIG(2)*SIG(5)-SIG(3)*SIG(4)) 
IF  (DET  .GT.  0.00000000001)  GOTO  30 
GOTO  20 
30   CONTINUE 

SIGINV(1)=(SIG(4)*SIG(6)-(SIG(5)**2))/DET 

SIGINV(2)=(SIG(3)*SIG(5)-SIG(2)*SIG(6))/DET 

SIGINV(3)=(SIG(2)*SIG(5)-SIG(3)*SIG(4))/DET 

SIGINV(4)=(SIG(1)*SIG(6)-(SIG(3)**2))/DET 

SIGINV(5)=(SIG(3)*SIG(2)-SIG(1)*SIG(5))/DET 

SIGINV(6)=(SIG(1)*SIG(4)-(SIG(2)**2))/DET 

DENOM=(SQRT(2.0*PI)**3)*SQRT(DET) 

CALL  FINDSEP 

RETURN 

END 


SUBROUTINE  MXTMLT(M1 ,M2, ANS, I , J,K) 

*  PERFORMS  MATRIX  MULTIPLICATION  M1 ( I , J)*M2( J,K)  TO  PRODUCE 

*  OUTPUT  MATRIX  ANS(I,K) 

INTEGER  L,ROW,COLfI,J,K 
REAL  M1(I,J)#M2(J,K),ANS(I,K) 
DO  10  ROW=1,I 
DO  20  COL=1,K 

ANS(ROW,COL)=0.0 
20     CONTINUE 
10  CONTINUE 

DO  30  ROW=1,I 
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DO  40  COL=1,K 
DO  50  L=1,J 

ANS(ROW,COL)=ANS(ROW,COL)+M1(ROW,L)*M2(L,COL) 
50        CONTINUE 
40     CONTINUE 
30  CONTINUE 
RETURN 
END 


SUBROUTINE  SEPEST(DT1 ,DT2,DT3,SEPM) 

PRODUCES  THE  9  BY  9  MATRIX  OF  PARTIALS  OF  SEP 

INCLUDE  'COM  DEF' 

REAL  DT1(10,9),DT2(22,10),DT3(9,22),SEPM(9,9),TEMP(22,9) 

CALL  MXTMLT(DT2,DT1,TEMP,22,10,9) 

CALL  MXTMLT(DT3,TEMP,SEPM,9,22,9) 

RETURN 

END 

SUBROUTINE  SETDAT 

SETS  THE  INPUT  VALUES  OF  MU,  SIG  AND  SEP  TO  BE  HELD  CONSTANT 
THROUGHOUT  THE  LOOPS  OF  THE  PROGRAM. 

INCLUDE  'COM  DEF' 
INCLUDE  'PMUCOM  DEF' 
INTEGER  I, J 
DO10  1=1,3 

MUSET(I)=MU(I) 
10  CONTINUE 
SEPSET=R 
RETURN 
END 


SUBROUTINE  TRIALS(N,MPT,TMU,TSIG) 

*  DRAWS  N  RANDOM  SAMPLES  FROM  NORMAL (MU, SIGMA)  AND  DETERMINES 

*  ESTIMATES  FOR  MU  AND  SIGMA  BASED  ON  THESE  SAMPLES. 

INCLUDE  'COM  DEF' 
INCLUDE  'PMUCOM  DEF' 
INTEGER  I,N,J 

REAL  MPT(3,3),RVAR(50,3),RSQ(50,3),RMLT(50,3),RMN(6),OP(3) 
C     ,NUMFAC,TSIG(6),TMU(3) 
EXTERNAL  DRNMVN,RNSET,DMACH,DRNUNF 
DO  10  1=1, N 

CALL  DRNMVN(1,3,MPT,3,0P,1) 

RVAR( I , 1 )=0P( 1 )+MUSET( 1 )/NUMSIM 

RVAR(I,2)=OP(2)+MUSET(2)/NUMSIM 

RVAR(I,3)=OP(3)+MUSET(3)/NUMSIM 

RSQ(I,1)=RVAR(I,1)**2 

RSQ(I,2)=RVAR(I,2)**2 

RSQ(I,3)=RVAR(I,3)**2 

RMLT( I , 1 )=RVAR( I , 1 )*RVAR( 1,2) 

RMLT( I ,2)=RVAR( I , 1 )*RVAR( 1,3) 

RMLT( I ,3)=RVAR( I ,2)*RVAR( 1,3) 
10  CONTINUE 

TMU(1)=FMN(RVAR,1,N) 

TMU(2)=FMN(RVAR,2,N) 

TMU(3)=FMN(RVAR,3,N) 

RMN(1)=FMN(RSQ,1,N) 

RMN(2)=FMN(RSQ,2,N) 

RMN(3)=FMN(RSO,3,N) 

RMN(4)=FMN(RMLT,1,N) 
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RMN(5)=FMN(RMLT,2,N) 

RMN(6)=FMN(RMLT,3,N) 

NUMFAC=1.0 

TSI G( 1 )=(RMN( 1 ) - TMU( 1 )**2)*NUMFAC 

TSIG(2)=(RMN(4)-TMU(1)*TMU(2)>*NUMFAC 

TSIG(3)=(RMN(5)-TMU(1)*TMU(3))*NUMFAC 

TSIG(4)=(RMN(2)-TMU(2)**2)*NUMFAC 

TSIG(5)=(RMN<6)-TMU(2)*TMU(3))*NUMFAC 

TSIG(6>=(RMN(3)-TMU(3)**2)*NUMFAC 

RETURN 

END 


REAL  FUNCTION  FMN(VEC,COL,N) 

DETERMINES  THE  MEAN  OF  N  ITEMS  IN  COLUMN  COL  OF  VECTOR  VEC 

INTEGER  I.N.COL 
REAL  VEC(50,3),TOT 
TOT =0.0 
DO  10  1=1, N 

TOT=TOT+VEC(I,COL> 
10  CONTINUE 

FMN=TOT/REAL(N) 

RETURN 

END 


REAL  FUNCTION  G(X) 

REAL  X 

G=0.0 

RETURN 

END 


REAL  FUNCTION  H(X) 

REAL  X 

H=6. 2831854 

RETURN 

END 


SUBROUTINE  WVALS 


*  DETERMINES  VALUES  OF  W,  CC,  CS,  SC  AND  SS  FUNCTIONS  USING  THE 

*  IMSL  SUBROUTINE  QTUODQ. 

REAL  G,H,H1,H2,ERREST 

INCLUDE  'COM  DEF' 

INCLUDE  'WVEC  DEF' 

EXTERNAL  G,H,DTWODQ,CC111 1 , CC131 1 ,CC31 1 1 ,CS1 111 ,CS131 1 , CS3111 , 
C  SC1 100, SC1200,SC1300,SC2111,SC3121, 882111,853112,553121, 
C  CC1113,CC1131,CC1211,CC1331,CC1411,CC1511,CC1531,CC2111, 
C  CC3113,CC3131,CC5111,CC5113,CC5131,CS1113,CS1131,CS1211, 
C  CS1331,CS1411,CS1511,CS1531,CS2111,CS3113,CS3131,CS5111, 
C  CS5113,CS5131,SC1121,SC1212,SC1221,SC1321,SC1400,SC1412, 
C  SC1421,SC1500,SC1521,SC3100,SC4111,SC4113,SC4131,SC5112, 
C  SC5114,SC5121,SC5141,SS1112,SS1121,SS1212,SS1221,SS1312, 
C  SS1321,SS1412,SS1421,SS1512,SS1521,SS4111,SS4113,SS4131, 
C       SS5112,SS5114,SS5121,SS5141 

B=PI 

CALL  DTWODQ  (CC1 1 1 1 ,A,B,G,H,ERRABS,ERRREL, IRULE,W(1 ),ERREST) 

CALL  DTWODQ  (CC131 1 ,A,B,G,H,ERRABS,ERRREL, IRULE,W(2),ERREST) 
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CALL 

DTWOOQ 

(CS1111,A,B,G,H, 

ERRABS.ERRREL, 

IRULE, 

W(3),ERREST) 

CALL 

DTWCOO 

(CS1311,A,B,G,H, 

ERRABS.ERRREL, 

IRULE, 

W(4),ERREST) 

CALL 

DTWOOQ 

<SC1100,A,B,GfH, 

ERRABS.ERRREL, 

IRULE, 

W(5),ERREST) 

CALL 

DTWOOQ 

(SC1200,A,B,G,H, 

ERRABS.ERRREL, 

IRULE, 

W(6),ERREST) 

CALL 

DTWOOQ 

(SC1300,A,B,G,H, 

ERRABS.ERRREL, 

IRULE, 

W(7),ERREST) 

CALL 

DTWOOQ 

(SC2111,A,B,G,H, 

ERRABS.ERRREL, 

IRULE, 

W(8),ERREST) 

CALL 

DTWOOQ 

(SC3121,A,B,G,H, 

ERRABS.ERRREL, 

IRULE, 

W(9),ERREST) 

CALL 

DTWOOQ 

(SS2111,A,B,G,H, 

ERRABS.ERRREL, 

IRULE, 

W(10),ERREST) 

CALL 

DTWOOQ 

(SS3112,A,B,G,H, 

ERRABS.ERRREL, 

IRULE, 

W(11),ERREST) 

CALL 

DTWOOQ 

(SS3121,A,B,G,H, 

ERRABS.ERRREL, 

IRULE, 

W(12),ERREST) 

CALL 

DTWOOQ 

(CC1113,A,B,G,H, 

ERRABS.ERRREL, 

IRULE, 

CC(D.ERREST) 

CALL 

DTWOOQ 

(CC1131,A,B,G,H, 

ERRABS.ERRREL, 

IRULE, 

CC(2),ERREST) 

CALL 

DTWOOQ 

(CC1211,A,B,G,H( 

ERRABS.ERRREL, 

IRULE, 

CC(3),ERREST) 

CALL 

DTWOOQ 

<CC1331,A,B,G,H, 

ERRABS.ERRREL, 

IRULE, 

CC(5),ERREST) 

CALL 

DTWOOQ 

(CC1411,A,B,G,H, 

ERRABS.ERRREL, 

IRULE, 

CC(6),ERREST) 

CALL 

DTWOOQ 

(CC1511,A,B,G,H, 

ERRABS.ERRREL, 

IRULE, 

CC(8),ERREST) 

CALL 

DTWOOQ 

(CC1531,A,B,G,H, 

ERRABS.ERRREL, 

IRULE, 

CC(9),ERREST) 

CALL 

DTWOOQ 

<CC2111,A,B,G,H, 

ERRABS.ERRREL, 

IRULE, 

CC(IO).ERREST) 

CALL 

DTWOOQ 

<CC3111,A,B,G,H, 

ERRABS.ERRREL, 

IRULE, 

CC(17),ERREST) 

CALL 

DTWOOQ 

(CC3113,A,B,G,H( 

ERRABS.ERRREL 

IRULE, 

CC(11),ERREST) 

CALL 

DTWOOQ 

(CCSISI.A.B.G.H, 

ERRABS.ERRREL 

IRULE 

CC(12),ERREST) 

CALL 

DTWODQ 

<CC5111,A,B,G,H 

ERRABS.ERRREL 

IRULE 

CC(14),ERREST) 

CALL 

DTWOOQ 

(CC5113,A,B,G,H 

ERRABS.ERRREL 

IRULE 

CC(15),ERREST) 

CALL 

DTWODQ 

(CC5131,A,B,G,H 

ERRABS.ERRREL 

IRULE 

CC(16),ERREST) 

CALL 

DTWOOQ 

(CS1113,A,B,G,H 

ERRABS.ERRREL 

IRULE 

CS(D.ERREST) 

CALL 

DTWODQ 

<CS1131,A,B,G,H 

ERRABS.ERRREL 

IRULE 

CS(2),ERREST) 

CALL 

DTWODQ 

<CS1211,A,B,G,H 

ERRABS.ERRREL 

IRULE 

CS(3),ERREST) 

CALL 

DTWODQ 

(CS1331,A,B,G,H 

ERRABS.ERRREL 

IRULE 

CS(5),ERREST) 

CALL 

DTWOOQ 

(CS1411,A,B,G,H 

ERRABS.ERRREL 

IRULE 

,CS(6),ERREST) 

CALL 

DTWODQ 

(CS1511,A,B,G,H 

ERRABS.ERRREL 

, IRULE 

,CS(8),ERREST) 

CALL 

DTWODQ 

(CS1531,A,B,G,H 

ERRABS.ERRREL 

, IRULE 

,CS(9),ERREST) 

CALL 

DTWODQ 

<CS2111,A,B,G,H 

.ERRABS.ERRREL 

, IRULE 

,CS(10),ERREST) 

CALL 

DTWOOQ 

(CS3111,AfB,G,H 

.ERRABS.ERRREL 

, IRULE 

,CS(17),ERREST) 

CALL 

DTWODQ 

<CS3113,A,B,G,H 

.ERRABS.ERRREL 

, IRULE 

,CS(18),ERREST) 

CALL 

DTWOOQ 

(CS3131,A,B,G,H 

.ERRABS.ERRREL 

, IRULE 

,CS(11),ERREST) 

CALL 

DTWODQ 

<CS5111,A,B,G,H 

.ERRABS.ERRREL 

, IRULE 

,CS(13),ERREST) 

CALL 

DTWODQ 

(CS5113,A,B,G,H 

.ERRABS.ERRREL 

, I  RULE 

,CS(14),ERREST) 

CALL 

DTWOOQ 

(CS5131,A,B,G,H 

.ERRABS.ERRREL 

, IRULE 

,CS(16),ERREST) 

CALL 

DTWODQ 

(SC1121,A,B,G,H 

.ERRABS.ERRREL 

, IRULE 

,SC(1),ERREST) 

CALL 

DTWODQ 

(SC1212,A,B,G,H 

.ERRABS.ERRREL 

, IRULE 

,SC(18),ERREST) 

CALL 

DTWODQ 

<SC1221,A,B,G,H 

.ERRABS.ERRREL 

, IRULE 

,SC(2),ERREST) 

CALL 

DTWODQ 

(SC1321,A,B,G,H 

.ERRABS.ERRREL 

, IRULE 

,SC(3),ERREST) 

CALL 

DTWODQ 

(SCKOO.A.B.G.H 

.ERRABS.ERRREL 

, IRULE 

,SC(4),ERREST) 

CALL 

DTWODQ 

(SC1412,A,B,G,H 

.ERRABS.ERRREL 

, IRULE 

,SC(19),ERREST) 

CALL 

DTWODQ 

(SC1421,A,B,G,H 

.ERRABS.ERRREL 

, IRULE 

,SC(5),ERREST) 

CALL 

DTWOOQ 

(SC1500,A,B,G,H 

.ERRABS.ERRREL 

, IRULE 

,SC(6),ERREST) 

CALL 

DTWODQ 

(SC1521,A,B,G,H 

.ERRABS.ERRREL 

, IRULE 

,SC(7),ERREST) 

CALL 

DTWOOQ 

(SC3100,A,B,G,H 

.ERRABS.ERRREL 

.IRULE 

,SC(9),ERREST) 

CALL 

DTWODQ 

(SC4111,A,B,G,H 

.ERRABS.ERRREL 

, IRULE 

,SC(10),ERREST) 

CALL 

DTWODQ 

(SC4113,A,B,G,H 

.ERRABS.ERRREL 

, IRULE 

,SC(11),ERREST) 

CALL 

DTWODQ 

(SC4131,A,B,G,H 

.ERRABS.ERRREL 

, IRULE 

,SC(12),ERREST) 

CALL 

DTWODQ 

(SC5112,A,BfG,H 

.ERRABS.ERRREL 

, IRULE 

,SC(13),ERREST) 

CALL 

DTWODQ 

(SC5114,A,B,G,H 

.ERRABS.ERRREL 

, IRULE 

,SC(14),ERREST) 

CALL 

DTWODQ 

(SC5121,A,B,G,H 

.ERRABS.ERRREL 

, IRULE 

,SC(15),ERREST) 

CALL 

DTWODQ 

(SC5141,A,B,G,H 

.ERRABS.ERRREL 

, I  RULE 

,SC(16),ERREST) 

CALL 

DTWODQ 

(SS1112,A,B,G,H 

.ERRABS.ERRREL 

, I  RULE 

.SS(D.ERREST) 

CALL 

DTWODQ 

(SS1121,A,B,G,H 

.ERRABS.ERRREL 

, I  RULE 

,SS(2),ERREST) 

CALL 

DTWODQ 

(SS1212,A,B,G,H 

.ERRABS.ERRREL 

, I  RULE 

,SS(3),ERREST) 

CALL 

DTWODQ 

<SS1221,A,B,G,H 

.ERRABS.ERRREL 

, IRULE 

,SS(4),ERREST) 

CALL 

DTWODQ 

<SS1312,A,B,G,H 

.ERRABS.ERRREL 

, IRULE 

,SS(5),ERREST) 

CALL 

DTWODQ 

<SS1321,A,B,G,H 

.ERRABS.ERRREL 

, IRULE 

,SS(6),ERREST) 

CALL 

DTWODQ 

(SS1412,A,B,G,H 

.ERRABS.ERRREL 

, IRULE 

,SS(7),ERREST) 

CALL 

DTWODQ 

(SS1421,A,B,G,H 

.ERRABS.ERRREL 

, IRULE 

,SS(8),ERREST) 

CALL 

DTWODQ 

<SS1512,A,B,G,H 

.ERRABS.ERRREL 

.IRULE 

,SS(9),ERREST) 

CALL 

DTWODQ 

(SS1521,A,B,G,H 

.ERRABS.ERRREL 

, IRULE 

,SS(10),ERREST) 

CALL 

DTWODQ 

(SS4111,A,B,G,H 

.ERRABS.ERRREL 

, IRULE 

,SS(12),ERREST) 

CALL 

DTWODQ 

(SS4113,A,B,G,H 

.ERRABS.ERRREL 

, IRULE 

,SS(13),ERREST) 

CALL 

DTWODQ 

(SS4131,A,B,G,H 

.ERRABS.ERRREL 

, IRULE 

,SS(14),ERREST) 
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CALL  DTWCOQ  (SS5112,A,B,G,H,ERRABS,ERRREL, IRULE,SS(15),ERREST) 
CALL  DTWOOQ  (SS51K, A, B,G,H,ERRABS, ERRREL, IRULE,SS(16),ERREST) 
CALL  DTWOOQ  (SS5121 , A,B,G,H,ERRABS, ERRREL, IRULE,SS(17),ERREST) 
CALL  DTWOOQ  (SS5K1 , A, B,G,H,ERRABS, ERRREL, IRULE, SS(18),ERREST) 
RETURN 
END 

SUBROUTINE  FINDSEP 

BASED  ON  PL- 1  PROGRAM  RAP,  PRODUCES  SEP  BASED  ON  MU  AND  SIGMA 
FINDSEP  DIAGONALIZES  SIGMA  BASED  ON  EIGENVALUES  AND  THEN  CALLS 
SUBROUTINE  ITERATION,  WHICH  DETERMINES  SEP 

INCLUDE  'COM  DEF' 

REAL  EIGVLS(3),EIGMAT(3,3),Q1(3,3),NMEAN(3),NEWSIG(3,3) 
C    ,TEMP(3,3),EIGNML(3,3) 

INTEGER  I, J 

EXTERNAL  DEVCSF 

CALL  DEVCSF(3,MSIG,3,EIGVLS,EIGMAT,3) 

CALL  ORTHO(EIGMAT,EIGNML) 

DO  11  1=1,3 
DO  21  J=1,3 

Q1(I,J)=EIGNML(J,I) 
21     CONTINUE 
11   CONTINUE 

CALL  MXTMLT(Q1,MU,NMEAN,3,3,1) 

CALL  MXTMLT(MSIG,EIGNML,TEMP,3,3,3) 

CALL  MXTMLT(Q1,TEMP,NEWSIG,3,3,3) 

CALL  ITERATION(NMEAN.NEWSIG) 

RETURN 

END 


SUBROUTINE  ORTHO  (A1.A2) 
ORTHONORMALIZES  MATRIX  OF  EIGENVECTORS  OF  SIGMA 

REAL  A1(3,3),A2(3,3),L,TEMP,TEMP2 

INTEGER  I 

L=SQRT((A1 (1 , 1 )**2)+(A1 (2,1 )**2)+<A1 (3, 1 )**2) ) 

DO  10  1=1,3 

A2(I,1)=A1(I,1)/L 
10  CONTINUE 

TEMP=A2( 1 , 1 )*A1 ( 1 ,2)+A2<2, 1 )*A1 (2,2)+A2(3, 1 )*A1 (3,2) 

DO  20  1=1,3 

A2(I,2)=A1(I,2)-TEMP*A2(I,1) 
20  CONTINUE 

L=SQRT(A2(1,2)**2+A2(2,2)**2+A2(3,2)**2) 

DO  30  1=1,3 

A2(I,2)=A2(I,2)/L 
30  CONTINUE 

TEMP=A2(1,1)*A1(1,3)+A2(2,1)*A1(2,3)+A2(3,1)*A1(3,3) 

TEMP2=A2(1,2)*A1(1,3)+A2(2,2)*A1(2,3)+A2(3,2)*A1(3,3) 

DO  40  1=1,3 

A2(I,3)=A1(I,3)-TEMP*A2(I,1)-TEMP2*A2(I,2) 
40  CONTINUE 

L=SQRT(A2(1,3)**2+A2(2,3)**2+A2(3,3)**2) 

DO  50  1=1,3 

A2(I,3)=A2(I,3)/L 
50  CONTINUE 

RETURN 

END 
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SUBROUTINE  ITERATION(MEAN, SIGMA) 

GENERATES  SEP  BASED  ON  MU  AND  SIGMA  FOR  DIAGONALIZED  SIGMA 

INCLUDE  'COM  DEF' 
INCLUDE  'SEPCOM  DEF' 

REAL  TOL,MEAN(3),SIGMA(3,3),MOM1,TRACE(3),SQSIG(3,3),CUSIG<3,3) 
C      ,TEMP(3),C2,DOF,RLOW,TPLOW,RHIGH,TPHIGH,CRATIO, INTEGRAL 
C      ,MOM2,MOM3 
TOL=0.001 

CXX=1.0/SIGMA(1,1) 
CYY=1.0/SIGMA(2,2) 
CZZ=1.0/SIGMA(3,3) 
CSANT=QUPITW*SQRT(CXX*CYY*CZZ) 

C=(CXX*(MEAN(1)**2)+CYY*(MEAN(2)**2)+CZZ*(MEAN(3)**2))/2.0 
CALL  MXTMLT(MEAN, MEAN, MOM1, 1,3,1) 
CALL  MXTMLT(SIGMA,SIGMA,SQSIG,3,3,3) 
CALL  MXTMLT(SIGMA,MEAN,TEMP,3,3,1) 
CALL  MXTMLT(MEAN, TEMP, MOM2, 1,3,1) 
CALL  MXTMLT(SIGMA,SQSIG,CUSIG,3,3,3) 
CALL  MXTMLT(SQSIG,MEAN,TEMP,3,3,1) 
CALL  MXTMLT (MEAN, TEMP, MOM3, 1,3,1) 
DO  10  1=1,3 

TRACE(I)=0.0 

MN(I)=MEAN(I) 
10  CONTINUE 
DO  20  1=1,3 

TRACE( 1 )=TRACE( 1 )+SIGMA( 1,1) 

TRACE(2)=TRACE(2)+SQSIG( 1,1) 

TRACE(3)=TRACE(3)+CUSIG( 1,1) 
20  CONTINUE 

MOM1=MOM1+TRACE(1) 
MOM2=2.0*(TRACE(2)+2.0*MOM2) 
MOM3=8.0*(TRACE(3)+3.0*MOM3) 
BETA=(MOM3**2)/(MOM2**3) 
DOF=8.0/BETA 

C2=DOF*(1.0-(2.0/(9.0*DOF)))**3 

RADIUS=SQRT(ABS(SQRT(ABS(MOM2/(2.0*DOF)))*(C2-DOF)+MOM1)) 
RLOU=0.95*RADIUS 
TPLOW=EVALT(RLOU) 
40  CONTINUE 

IFCTPLOW  .GT.  0.5)  THEN 

RLOW=0.9*RLOW 

TPLOW=EVALT(RLOW) 
GOTO  40 
END  IF 

RHIGH=1.105*RLOU 
TPHIGH=EVALT(RHIGH) 
50  CONTINUE 

IF  (TPHIGH  .LT.  0.5)  THEN 

RHIGH=1.1*RHIGH 

TPHIGH=EVALT(RHIGH) 
GOTO  50 
END  IF 

CRATIO=(0.5-TPLOW)/(TPHIGH-TPLOW) 
RADIUS=RLOW+(RHIGH-RLOW)*CRATIO 
ICNT=0 
60  CONTINUE 
ICNT=ICNT+1 
IF((RHIGH-RLOW)/RHIGH  .GT.  TOL)  THEN 

INTEGRAL=EVALT(RADIUS) 

I F(ABS( INTEGRAL  -  0.5)  .LT.  TOL)  GOTO  75 

IFONTEGRAL  .GT.  0.5)  THEN 
RHIGH  =  RADIUS 
TPHIGH=INTEGRAL 

ELSE 

RLOU=RADIUS 
TPLOW=INTEGRAL 
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END  IF 

RADIUS=RLOW+(RHIGH-RLOW)*CRATIO 
IF  (ICNT  .GE.  150)  THEN 
PRINT*, 'COUNT  EXCEEDED1 
GOTO  75 
END  IF 
GOTO  60 
ENDIF 
75  CONTINUE 
R=RADIUS 
RETURN 
END 


REAL  FUNCTION  EVALT(RAD) 

EVALUATES  THE  AREA  COVERAGE  BY  TRIVARIATE  NORMAL  FOR  INPUT  R 

INCLUDE  'COM  DEF' 
INCLUDE  'SEPCOM  DEF' 
REAL  ANS, LOWER .UPPER, RAD 
EXTERNAL  RAPFC3,RAPLOW,RAPUP 
RADIUS=RAD 
LOWER=0.0 
UPPER=6. 2831853 

CALL  DTWODQ(RAPFC3, LOWER, UPPER, RAPLOW,RAPUP,ERRABS,ERRREL, IRULE, 
C  ANS.ERREST) 

EVALT=ANS 

RETURN 

END 


REAL  FUNCTION  RAPLOW(X) 

REAL  X 

RAPLOW=0.0 

RETURN 

END 


REAL  FUNCTION  RAPUP(X) 

REAL  X 

RAPUP=3. 1415927 

RETURN 

END 


REAL  FUNCTION  RAPFC3(THETA,PHI ) 

FUNCTION  USED  TO  FIND  SEP 

REAL  CP,SP,CT,ST,AA,SA,SAM1,BB,B2DA,ERFARG,EXPARG,THETA,PHI 

INCLUDE  'COM  DEF' 

INCLUDE  'SEPCOM  DEF' 

EXTERNAL  ERF 

CP=COS(PHI) 

CT=COS(THETA) 

SP=SIN(PHI) 

ST=SIN(THETA) 

AA=(CXX*SP*SP*CT*CT+CYY*SP*SP*ST*ST+CZZ*CP*CP)/2.0 

SA=SQRT(AA) 

SAM1=1.0/SA 

BB=(CXX*MN(1)*SP*CT+CYY*MN(2)*SP*ST+CZZ*MN(3)*CP)/2.0 

B2DA=BB*BB/AA 

ER  FARG=SA*RAD I US+SAM1 *BB 

EXPARG= - AA*RAD I US*RAD I  US  -  2 . 0*BB*RAD I  US  - B2DA 

IF  (EXPARG  .GT.  -32.0)  THEN 

F=ERF(ERFARG)*SAM1*(0.5+B2DA)-EXP(EXPARG)*(RADIUS-BB/AA) 
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C        *SQPIINV 
ELSE 

F=ERF(ERFARG)*SAM1*(0.5+B2DA) 
ENDIF 

ERFARG=SAM1*BB 
EXPARG=-B2DA 
IF  (EXPARG  .GT.  -32.0)  THEN 

F=F-ERF(ERFARG)*SAM1*(0.5+B2DA)-EXP(EXPARG)*(BB/AA)*SQPIINV 
ELSE 

F=F-ERF(ERFARG)*SAM1*(0.5+B2DA) 
ENDIF 

EXPARG=B2DA-C 
IF  (EXPARG  .GT.  -32.0)  THEN 

F=F*CSANT*EXP(EXPARG)*SP/AA 

IF(F  .LT.  0.000000000000001)  F=0.0 
ELSE 

F=0.0 
ENDIF 
RAPFC3=F 
RETURN 
END 


SUBROUTINE  DT1VAL  (DT1 , SEPMU,SEPSIG) 

SETS  THE  VALUES  OF  THE  MATRIX  DT1  AND  THE  VECTORS  SEPMU  AND 
SEPSIG. 

INCLUDE  'COM  DEF' 
INCLUDE  'WVEC  DEF1 
REAL  DT1(10,9),SEPMU(3),SEPSIG(6) 
INTEGER  I, J 
DO  100  1=1,9 
DO  200  J=1,9 

IF  (I  .NE.  J)  THEN 

DT1(I,J)=0.0 
ELSE 

DT1(I,J)=1.0 
ENDIF 
200    CONTINUE 
100  CONTINUE 

DT1(10,1)=U(8)/W(5) 
DT1(10,2)=W(10)/W(5) 
DT1(10,3)=W(6)/(2.0*W(5)) 

DT1(10,4)=(R/<8.0*W(5)))*(4.0*SIGINV<1)*W(9)+2.0*SIGINV(2)*W(11) 
C         +SIGINV(3)*(W(1)-W(2)))-(W(8)/(2.0*W(5)))* 
C  (SIGINV(1)*MU(1)+SIGINV(2)*MU(2)+SIGINV(3)*MU(3)) 

DT1(10,5)=(R/(4.0*W(5)))*(2.0*SIGINV(1)*W(11)+4.0*SIGINV(2)*W(12) 
C         +SIGINV(3)*(W(3)-W(4)))-(W(10)/W(5))* 
C  <SIGINV(1)*MU(1)+SIGINV(2)*MU(2)+SIGINV(3)*MU(3)) 

DT1(10,6)=(R/(4.0*W(5)))*(SIGINV<1)*(W(1)-W(2))+SIGINV(2)*(W(3) 
C         -U(4))+SIGINV(3)*(U(7)+U(5)))-(W(6)/(2.0*W(5)))* 
C  (SIGINV(1)*MU(1)+SIGINV(2)*MU(2)+SIGINV(3)*MU(3)) 

DT1(10,7)=(R/(8.0*W(5)))*(2.0*SIGINV(2)*W(11)+4.0*SIGINV(4) 
C  *W(12)+SIGINV(5)*(W(3)-W(4)))-(W(10)/(2.0*W(5)))* 
C  (SIGINV(2)*MU(1)+SIGINV(4)*MU(2)+SIGINV(5)*MU(3)) 

DT1(10,8)=(R/(4.0*W(5)))*(SIGINV(2)*(W(1)-W(2))+SIGINV(4)*(W(3) 
C         -W(4))+SIGINV(5)*(W(7)+W(5)))-(W(6)/(2.0*W(5)))* 
C  (SIGINV(2)*MU(1)+SIGINV(4)*MU(2)+SIGINV(5)*MU(3)) 

DT1(10,9)=(R/(8.0*W(5)))*(SIGINV(3)*(W(1)-W(2))+SIGINV(5)*(W(3) 
C  -W(4))+SIG1NV(6)*(W(7)+W(5)))-(W(6)/(4.0*W(5)))* 

C  (SIGINV(3)*MU(1)+SIGINV(5)*MU(2)+SIGINV(6)*MU(3)) 

DO  10  1=1,3 

SEPMU(I)=DT1(10,I) 
10  CONTINUE 
DO  20  1=1,6 

SEPSIG(I)=DT1(10,I+3) 
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20  CONTINUE 
RETURN 
END 


SUBROUTINE  DT2VAL(DT2) 

*  SETS  THE  VALUES  OF  THE  MATRIX  DT2.   THE  VECTORS  A1  AND  B1  ARE 
'        IN  THIS  SUBROUTINE  AND  THEN  PASSED  TO  THE  SUBROUTINE  FIXMAT, 

*  WHICH  ACTUALLY  SETS  THE  VALUES  OF  DT2. 

INCLUDE  'COM  DEF' 

INCLUDE  'WVEC  DEF' 

REAL  DT2(22,10),SIGM(3,3),TEMP(3),ANS<3),B1<3),A1(6) 

INTEGER  I, J 

DO  100  1=1,10 

DO  200  J=1,10 

IF  (I  .ME.  J)  THEN 

DT2(I,J)=0.0 
ELSE 

DT2(I,J)=1.0 
END  IF 
200     CONTINUE 
100  CONTINUE 

SIGM(1/1)=SIGINV(1) 
SIGM(1,2)=SIGINV(2) 
SIGM(1,3)=SIGINV(3) 
SIGM(2,1)=SIGINV(2) 
SIGM(2,2)=SIGINV(4) 
SIGM(2,3)=SIGINV(5) 
SIGM(3,1)=SIGINV(3) 
SIGM(3,2)=SIGINV(5) 
SIGM(3,3)=SIGINV(6) 

B1(1)=SC(2)/2.0 

B1(2)=SS(3)/4.0 

B1(3)=CC(10) 

A1(1)=CC(2)-CC<12) 

A1(2)=2.0*(W(3)-CS(17)-CS(2)+CS(11)) 

A1(3)=2.0*(SC(1)-W(9)) 

A1(4)=U(1)-CC<17)-CC(2)+CC(12) 

A1(5)=SS(1)-W(11) 

A1(6)=CC(17) 
CALL  FIXMAT(A1,B1,11,DT2,SIGM) 

B1(1)=(SC(5)-SC(2))/2.0 

B1(2)=(SS(7)-SS(3))/4.0 

B1(3)=(CC(6)+CC(3))/2.0 

A1(1)=(2.0*CC(5)-CC(9)-CC(2))/4.0 

A1(2)=<2.0*W(4)-CS(8)-W(3)-2.0*CS(5)+CS(9)+CS(2))/2.0 

A1(3)=(SC(7)-SC(1))/2.0 

A1(4)=(2.0*U(2)-CC(8)-U(1)-2.0*CC(5)+CC(9)+CC(2))/4.0 

A1(5)=(SS(9)-SS(1))/4.0 

A1(6)=(CC(8)+2.0*W(2)+W(1))/4.0 
CALL  FIXMAT(A1,B1,12,DT2,SIGM) 

B1(1)=SS(3)/4.0 

B1(2)=SS(4)/2.0 

B1(3)=CS(10) 

A1(1)=W(3)-CS(17)-CS(2)+CS(11) 

A1(2)=2.0*(W(1)-CC(2)-CC(17)+CC(12)) 

A1(3)=SS(1)-W(11) 

A1(4)=CS(2)-CS(11) 

A1(5)=2.0*(SS(2)-U(12)) 

A1(6)=CS(17) 
CALL  FIXMAT(A1,B1,13,DT2,SIGM) 

B1(1)=(SS(7)-SS(3))/4.0 

B1(2)=(SS(8)-SS(4))/2.0 

B1(3)=(CS(6)+CS(3))/2.0 

A1(1)=(2.0*W(4)-2.0*CS(5)-CS(8)+CS(9)-W(3)+CS(2))/4.0 

A1(2)=(2.0*W(2)-2.0*CC(5)-CC(8)+CC(9)-U(1)+CC(2))/2.0 
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A1(3)=(SS(9)-SS(1))/4.0 

A1(4)=(2.0*CS(5)-CS(9)-CS(2))/4.0 

A1(5)=(SS(10)-SS(2>)/2.0 

Al(6)=(CS(8)+2.0*U(4)+U(3))/4.0 
CALL  FIXMAT(A1,B1,14,DT2,SIGM) 

B1(1)=W(8) 

B1(2)=W(10) 

B1(3)=W(6)/2.0 

A1(1)=W(9) 

A1(2)=W(11) 

A1(3)=2.0*(U(1)-CC(17>> 

A1(4)=W(12) 

A1(5)=2.0*(U(3)-CS(17)) 

A1(6)=W(5)-SC<9) 
CALL  FIXMAT(A1,B1,15,DT2,SIGM) 

B1(1)=(W(1)-W(2))/2.0 

B1(2)=(W(3)-W(4))/2.0 

B1(3)=(W(7)+W(5))/2.0 

A1(1)=(2.0*SC(2)-SC(5))/4.0 

A1(2)=(2.0*SS(3)-SS(7))/4.0 

A1(3)=(2.0*W(8)+CC(3)CC(6))/2.0 

A1(4)=(2.0*SS(4)-SS(8))/4.0 

A1(5)=(2.0*U(10)+CS(3)-CS(6))/2.0 

A1(6)=(2.0*W(6)+SC(4))/4.0 
CALL  FIXMAT(A1,B1,16,DT2,SIGM) 

B1(1)=(CC(3)-CC(6))/2.0 

B1(2)=(CS(3)CS(6))/2.0 

B1(3)=(SC(4)+W(6))/2.0 

A1(1)=(2.0*SC(3)-SC(1)-SC(7))/4.0 

A1(2)=(2.0*SS(5)-SS(1)-SS(9))/4.0 

A1(3)=(W(1)-CC(8)>/2.0 

A1(4)=(2.0*SS(6)-SS(2)-SS(10))/4.0 

A1(5)=(W(3)-CS(8))/2.0 

A1(6)=(SC(6)+2.0*W(7)+W(5))/4.0 
CALL  FIXMAT(A1,B1,17fDT2,SIGM) 

B1(D=W(9) 

B1(2)=W(11)/2.0 

B1<3)=U(1)-CC(17) 

A1(1)=SC(12) 

A1(2)=2.0*(SS(12)-SS(14)) 

A1(3)=(2.0*SC(2)-SC(5))/4.0 

A1(4)=SC(10)-SC(12) 

A1(5)=(2.0*SS(3)-SS(7))/8.0 

A1(6)=W(8)-SC(10) 
CALL  FIXMAT(A1,B1,18,DT2,SIGM) 

B1(1)=SC(12) 

B1(2)=SS(12)-SS(14) 

B1(3)=(2.0*SC(2)-SC(5))/8.0 

A1(1)=SC(16) 

A1(2)=(SS(16)+2.0*SS(15))/4.0 

A1(3)=2.0*(CC(2)-2.0*CC(12)+CC(16)) 

A1(4)=SS(17)-SS(18) 

A1(5)=2.0*(W(3)-CS(2)-2.0*CS<17)+2.0*CS(11)+CS(13)-CS(16)) 

A1(6)=W(9)-SC(15) 
CALL  FIXMAT(A1,B1,19,DT2,SIGM) 

B1(1)=W(11)/2.0 

B1(2)=W(12) 

B1(3)=W(3)-CS(17) 

A1(1)=SS(12)-SS(14) 

A1(2)=2.0*(SC(10)-SC(12)) 

A1(3)=(2.0*SS(3)-SS(7))/8.0 

A1(4)=SS(14) 

A1(5)=(2.0*SS(4)-SS(8))/4.0 

A1(6)=W(10)-SS(12) 
CALL  FIXMAT(A1,B1,20,DT2,SIGM) 

B1(1)=(SS(13)+SS(12))/2.0 

B1(2)=(SC(10)-SC(11))/2.0 

B1(3)=(2.0*SS(3)-SS(7))/8.0 
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A1(1)=(SS(16)+2.0*SS(15))/4.0 

A1(2)=(2.0*SC(15)-SC(14)-SC(13))/2.0 

A1(3)=CS(1)+U(3)-2.0*CS(18)-2.0*CS(17)+CS(14)+CS(13) 

A1(4)=<2.0*SS(15)-SS<16))/4.0 

A1(5)=W(1)-CC(1)-2.0*CC(17)+2.0*CC(11)+CC(14)-CC(15) 

A1(6)=W(11)-SS(15) 
CALL  FIXMAT(A1,B1,21,DT2,SIGM) 

B1(1)=SC(10)-SC(12) 

B1(2)=SS(14) 

B1(3)=(2.0*SS(4)-SS(8))/8.0 

A1(1)=SS(17)-SS(18) 

A1(2)=(2.0*SS(15)-SS(16))/4.0 

A1(3)=2.0*(W(1)-CC(2)-2.0*CC(17)+2.0*CC(12)+CC(14)-CC(16)) 

A1(4)=SS(18) 

A1(5)=2.0*(CS(2>-2.0*CS(11)+CS(16)) 

A1(6)=W(12)-SS(17) 
CALL  FIXMAT(A1,B1,22,DT2,SIGM) 
RETURN 
END 


SUBROUTINE  FIXMAT(A1 ,B1 , N,DT2,SIGM) 

*  SETS  THE  VALUES  OF  DT2,  GIVEN  A1,  B1,  SIGMA  INVERSE  (SIGM) 

*  AND  THE  ROW  NUMBER  TO  BE  SET  AS  INPUTS. 

INCLUDE  'WVEC  DEF' 
INCLUDE  'COM  DEF' 

REAL  DT2(22,10),SIGM(3,3),A1<6),B1(3),A2(6),A3<6),A4(3),MSUM(6) 
C       ,G(6),ANS<3) 
INTEGER  N,I 
DO  100  1=1,3 

A4(I)=B1(I)*R-MU(I)*W(N-10) 
100  CONTINUE 

CALL  MXTMLT(SIGM,A4,ANS, 3,3,1) 

DT2(N,1)=ANS(1) 

DT2(N,2)=ANS(2) 

DT2(N,3)=ANS(3) 

A2(1)=B1(1)*MU(1) 

A2(2)=B1(2)*MU(1)+B1(1)*MU(2) 

A2(3)=B1(3)*MU(1)+B1(1)*MU<3) 

A2(4)=B1(2)*MU(2) 

A2(5)=B1(3)*MU(2)+B1(2)*MU(3) 

A2(6)=B1(3)*MU(3) 

DO  200  1=1,6 

MSUM(I)=R*A1(I)-A2(I) 
200  CONTINUE 

CALL  MXTMLT(SIGINV,MSUM,ANS,  1,6,1) 

DT2(N,10)=-ANS(1) 

A3(1)=W(N-10)*(MU(1)**2) 

A3(2)=2.0*W(N-10)*MU(1)*MU(2) 

A3(3)=2.0*W(N-10)*MU(1)*MU(3) 

A3(4)=W(N-10)*(MU(2)**2) 

A3(5)=2.0*W(N-10)*MU(2)*MU(3) 

A3(6)=W(N- 10)*(MU(3)**2) 

DO  300  1=1,6 

MSUM( I )=(R**2)*A1 ( I ) -2.0*R*A2( I )+A3( I ) 
300  CONTINUE 

G(1)=SIGINV(1)*SIGINV(1) 

G(2)=SIGINV(2)*SIGINV(1) 

G(3)=SIGINV(3)*SIGINV(1) 

G(4)=SIGINV(4)*SIGINV(1)-SIG(6)/DET 

G(5)=SIGINV(5)*SIGINV(1)+SIG(5)/DET 

G(6)=SIGINV(6)*SIGINV(1)-SIG(4)/DET 

CALL  MXTMLT(G,MSUM,ANS, 1,6,1) 

DT2(N,4)=(ANS(1)-W(N-10)*SIGINV(1))/2.0 

G(1)=SIGINV(1)*SIGINV(2) 

G(2)=SIGINV(2)*SIGINV(2)+SIG(6)/(2.0*DET) 
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G(3)=SIGINV(3)*SIGINV(2)-SIG(5)/(2.0*DET) 

G(4)=SIGINV(4)*SIGINV(2) 

G(5)=SIGINV(5)*SIGINV(2)-SIG(3)/(2.0*DET) 

G(6)=SIGINV(6)*SIGINV(2)+SIG(2)/DET 

CALL  MXTMLT(G,MSUM,ANS,  1,6,1) 

DT2(N,5)=ANS(1)-W(N-10)*SIGINV(2) 

G(1)=SIGINV(1)*SIGINV(3) 

G(2)=SIGINV(2)*SIGINV(3)-SIG(5)/(2.0*DET) 

G(3)=SIGINV(3)*SIGINV(3)+SIG<4)/(2.0*DET) 

G(4)=SIGINV(4)*SIGINV(3)+SIG(3)/DET 

G(5)=SIGINV(5)*SIGINV(3)-SIG(2)/(2.0*DET) 

G(6)=SIGINV(6)*SIGJNV(3) 

CALL  MXTMLT(G,MSUM,ANS, 1,6,1) 

DT2(N,6)=ANS(1)-W(N-10)*SIGINV(3) 

G(1)=SIGINV(1)*SIGINV(4)-SIG(6)/DET 

G(2)=SIGINV(2)*SIGINV(4) 

G(3)=SIGINV(3)*SIGINV(4)+SIG(3)/DET 

G(4)=SIGINV(4)*SIGINV(4) 

G(5)=SIGINV(5)*SIGINV(4) 

G(6)=SIGINV(6)*SIGINV(4)-SIG(1)/DET 

CALL  MXTMLT(G,MSUM,ANS, 1,6,1) 

DT2(N,7)=(ANS(1)-U(N-10)*SIGINV(4))/2.0 

G(1)=SIGINV(1)*SIGINV(5)+SIG(5)/DET 

G(2)=SIGINV(2)*SIGINV(5)-SIG(3)/(2.0*DET) 

G(3)=SIGINV(3)*SIGINV(5)-SIG(2)/(2.0*DET) 

G(4)=SIGINV(4)*SIGINV<5) 

G(5)=SIGINV(5)*SIGINV(5)+SIG(1)/(2.0*DET) 

G(6)=SIGINV(6)*SIGINV(5) 

CALL  MXTMLT(G,MSUM,ANS, 1,6,1) 

DT2(N,8)=ANS(1)-W(N-10)*SIGINV(5) 

G(1)=SIGINV(1)*SIGINV(6)-SIG(4)/DET 

G(2)=S1GINV(2)*SIGINV(6)+SIG(2)/DET 

G(3)=SIGINV(3)*SIGINV(6) 

G(4)=SIGINV(4)*SIGINV(6)-SIG(1)/DET 

G(5)=SIGINV(5)*SIGINV(6) 

G(6)=SIGINV(6)*SIGINV(6) 

CALL  MXTMLT(G,MSUM,ANS, 1,6,1) 

DT2(N,9)=(ANS(1)-W(N-10)*SIGINV(6))/2.0 

RETURN 

END 


SUBROUTINE  DT3VAL(DT3) 

SETS  THE  VALUES  OF  THE  MATRIX  DT3 


INCLUDE  'COM  DEF' 

INCLUDE  'WVEC  DEF' 

REAL  DT3(9,22),U41 ,U42,U43,U51 ,U52,U61 ,U62,U71 ,U72,U73,U81 ,U82, 
C     U91.U92.U93 

INTEGER  I, J 

DO  100  1=1,3 
DO  200  J=1,22 
DT3(I,J)=0.0 
200    CONTINUE 
100  CONTINUE 

DT3(1,15)=-U(8)/(U<5)**2) 

DT3(1,18)=1.0/W(5) 

DT3(2,15)=-W(10)/(W(5)**2) 

DT3(2,20)=1.0/W(5) 

DT3(3,15)=-0.5*W(6)/(W(5)**2) 

DT3(3,16)=0.5/W(5) 

U41=R*(2.0*SIGINV(1)*W(9)+SIGINV(2)*U(11)+0.5*SIGINV(3)*(U(1) 
C       -W(2)))/(4.0*W(5)) 

U42=W(8)/(2.0*W(5)) 

U43=SIGINV(1)*MU(1)+SIGINV(2)*MU(2)+SIGINV(3)*MU(3) 

U4=U41-U42*U43 
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DT3(4,1)=-SIGINV(1)*U42 

DT3(4,2)=-SIGINV(2)*U42 

DT3(4,3)=-SIGINV(3)*U42 

U51=R*(SIGINV(1)*W(11)+2.0*SIGINV(2)*W(12)+0.5*SIGINV(3)*(W(3) 
C       -W(4)))/(2.0*W(5>) 

U52=W(10)/W(5) 

U5=U51-U52*U43 

DT3(5,1)=-SIGINV(1)*U52 

DT3(5,2)=-SIGINV(2)*U52 

DT3(5,3)=-SIGINV(3)*U52 

U61=R*(SIGINV(1)*(W(1)-W(2))+SIGINV(2)*(W(3)-W(4))+SIGINV(3) 
C       *(W(7)+W(5)))/(4.0*W(5)) 

U62=W(6)/(2.0*W(5)) 

U6=U61-U62*U43 

DT3(6,1)=-SIGINV(1)*U62 

DT3(6f2)=-SIGINV(2)*U62 

DT3(6,3)=-SIGINV(3)*U62 

U71=R*(SIGINV(2)*U(11)+2.0*SIGINV(4)*U(12)+0.5*SIGINV(5) 
C       *(W(3)-W(4)))/(4.0*W(5)) 

U72=W(10)/(2.0*W(5)) 

U73=SIGINV(2)*MU(1)+SIGINV(4)*MU(2)+SIGINV(5)*MU<3> 

U7=U71-U72*U73 

DT3(7,1)=-SIGINV(2)*U72 

DT3(7,2)=-SIGINV(4)*U72 

DT3(7,3)=-SIGINV(5)*U72 

U81=R*(SIGINV(2)*(W(1)-W(2))+SIGINV(4)*(W(3)-W(4))+SIGINV(5) 
C       *(U(7)+W(5)))/(4.0*W(5)> 

U82=U(6)/(2.0*W(5)) 

U8=U81-U82*U73 

DT3(8,1)=-SIGINV(2)*U82 

DT3(8,2)=-SIGINV(4)*U82 

DT3(8,3)=-SIGINV(5)*U82 

U91=R*(SIGINV(3)*(W(1)-W(2))+SIGINV(5)*(W(3)-W(4))+SIGINV(6) 
C       *(W(7)+W(5)))/(8.0*W(5)) 

U92=W(6)/(4.0*W(5)) 

U93=SIGINV(3)*MU(1)+SIGINV(5)*MU(2)+SIGINV(6)*MU(3) 

U9=U91-U92*U93 

DT3(9,1)=-SIGINV(3)*U92 

DT3(9,2)=-SIGINV(5)*U92 

DT3(9,3)=-SIGINV(6)*U92 

TEMP1=((R*W(11)/2.0)-MU(2)*W(8))/(2.0*W(5)) 

TEMP2=(0.25*R*(W(1)-W(2))-MU(3)*W(8))/(2.0*W(5)) 

TEMP3=(R*W(9)-MU(1)*W(8))/(2.0*W(5)) 

DT3(4,4)=-SIGINV(1)*U4 

DT3(4,5)=-2.0*SIGINV(2)*U4-(SIG(6)*TEMP1-SIG(5)*TEMP2)/DET 

DT3(4,6)=-2.0*SIGINV(3)*U4-(SIG(4)*TEMP2-SIG(5)*TEMP1)/DET 

DT3(4,7)=-SIGINV(4)*U4-(SIG(3)*TEMP2-SIG(6)*TEMP3)/DET 

DT3(4,8)=-2.0*SIGINV(5)*U4-(2.0*SIG(5)*TEMP3-SIG(3)*TEMP1-SIG(2) 
C        *TEMP2)/DET 

DT3(4,9)=-SIGINV(6)*U4-(SIG(2)*TEMP1-SIG(4)*TEMP3)/DET 

TEMP1=(R*W(12)-MU(2)*W(10))/W(5) 

TEMP2=(0.25*R*(W(3)-W(4))-MU(3)*W(10))/W(5) 

TEMP3=(R*W(11)/2.0-MU(1)*W(10))/W(5) 

DT3(5,4)=-SIGINV(1)*U5 

DT3(5,5)=-2.0*SIGINV(2)*U5-(SIG(6)*TEMP1-SIG(5)*TEMP2)/DET 

DT3(5,6)=-2.0*SIGINV(3)*U5-(SIG(4)*TEMP2-SIG(5)*TEMP1)/DET 

DT3(5,7)=-SIGINV(4)*U5-(SIG(3)*TEMP2-SIG(6)*TEMP3)/DET 

DT3(5f8)=-2.0*SIGINV(5)*U5-(2.0*SIG(5)*TEMP3-SIG(3)*TEMP1-SIG(2) 

C        *TEMP2)/DET 
DT3(5,9)=-SIGINV(6)*U5-(SIG(2)*TEMP1-SIG(4)*TEMP3)/DET 
DT3(7,4)=-SIGINV(1)*U7-(SIG(5)*TEMP2-SIG(6)*TEMP1)/(2.0*DET) 
DT3(7,5)=-2.0*SIGINV(2)*U7-(SIG(6)*TEMP3-SIG(3)*TEMP2)/(2.0*DET) 
DT3(7,6)=-2.0*SIGINV(3)*U7-(2.0*SIG(3)*TEMP1-SIG(5)*TEMP3-SIG(2) 

C        *TEMP2)/(2.0*DET) 
DT3(7,7)=-SIGINV(4)*U7 

DT3(7,8)=-2.0*SIGINV(5)*U7-(SIG(1)*TEMP2-SIG(3)*TEMP3)/(2.0*DET) 
DT3(7,9)=-SIGINV(6)*U7-(SIG(2)*TEMP3-SIG(1)*TEMP1)/(2.0*DET) 
TEMPl=(R*(W(3)-W(4))/2.0-MU(2)*W(6))/(2.0*W(5)) 
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400 
300 


TEMP2=(R*(W(7)+U(5))/2.0-MU(3)*U(6))/(2.0*U(5)) 

TEMP3=(R*(W(1)-W(2))/2.0-MU(1)*W(6))/(2.0*W(5)) 

DT3(6,4)=-SIGINV(1)*U6 

DT3(6,5)=-2.0*SIGINV(2)*U6-(SIG(6)*TEMP1-SIG(5)*TEMP2)/DET 

DT3(6,6)=-2.0*SIGINV(3)*U6-(SIG(4)*TEMP2-SIG(5)*TEMP1)/DET 

DT3(6,7)=-SIGINV(4)*U6-(SIG(3)*TEMP2-SIG(6)*TEMP3)/DET 

DT3(6,8)=-2.0*SIGINV(5)*U6-(2.0*SIG(5)*TEMP3-SIG(3)*TEMP1-SIG(2) 

C        *TEMP2)/DET 

DT3(6,9)=-SIGINV(6)*U6-(SIG(2)*TEMP1-SIG(4)*TEMP3)/DET 
DT3(8,4>=-SIGINV(1)*U8-(SIG(5)*TEMP2-SIG(6)*TEMP1)/DET 
DT3(8,5)=-2.0*SIGINV(2)*U8-(SIG(6)*TEMP3-SIG(3)*TEMP2)/DET 
DT3(8,6)=-2.0*S1GINV(3)*U8-(2.0*SIG(3)*TEMP1-S1G(5)*TEMP3-SIG(2) 

C        *TEMP2)/DET 
DT3(8,7)=-SIGINV(4)*U8 

DT3(8,8)=-2.0*SIGINV(5)*U8-(SIG(1)*TEMP2-SIG(3)*TEMP3)/DET 
DT3(8,9)=-SIGINV(6)*U8-(SIG(2)*TEMP3-SIG(1)*TEMP1)/DET 
DT3(9,4)=-SIGINV(1)*U9-(SIG(5)*TEMP1-SIG(4)*TEMP2)/(2.0*DET) 
DT3(9,5)=-2.0*SIGINV(2)*U9-(2.0*SIG(2)*TEMP2-SIG(3)*TEMP1-SIG(5) 

C        *TEMP3)/(2.0*DET) 

DT3(9,6)=-2.0*SIGINV(3)*U9-(SIG(4)*TEMP3-SIG(2)*TEMP1)/(2.0*DET) 
DT3(9,7)=-SIGINV(4)*U9-(SIG(3)*TEMP3-SIG(1)*TEMP2)/(2.0*DET) 
DT3(9,8)=-2.0*SIGINV(5)*U9-(SIG(1)*TEMP1-SIG(2)*TEMP3)/(2.0*DET) 
DT3(9,9)=-SIGINV(6)*U9 
DT3(4,10)=U41/R 
DT3(5,10)=U51/R 
DT3(6,10)=U61/R 
DT3(7,10)=U71/R 
DT3(8,10)=U81/R 
DT3(9,10)=U91/R 
DO  300  1=4,9 

DO  400  J=11,22 

DT3(I,J>=0.0 
CONTINUE 
CONTINUE 

DT3(4,11)=(R*SIGINV(3))/(8.0*W(5)) 
,12)=-DT3(4,11) 


DT3(4 
DT3(4 
DT3(4 
DT3(4 
DT3(4 
DT3(5 
DT3(5 
DT3(5 
DT3(5 
DT3(5 
DT3(5 
DT3(6 
DT3(6 
DT3(6 
DT3(6 
DT3(6 
DT3(6 
DT3(6 
DT3(7 
DT3(7 
DT3(7 
DT3(7 
DT3(7 
DT3(7 
DT3(8 
DT3(8 
DT3(8 
DT3(8 
DT3(8 
DT3(8 
DT3(8 
DT3(9 
DT3(9 


15)=-U4/W(5) 

18)=-U43/(2.0*W(5)) 

19)=(R*SIGINV(1))/(2.0*W(5)) 

21)=(R*SIGINV(2))/(4.0*W(5)) 

13)=(R*SIGINV(3))/(4.0*W(5)) 

14)=-DT3(5,13) 

15)=-U5/W(5) 

20)=-U43/W(5) 

21)=DT3(4,19) 

22)=(R*SIGINV(2))/W(5) 

11)=(R*SIGINV(1))/(4.0*W(5)) 

12)=-DT3(6,11) 

13)=DT3(4,21) 

14)=-DT3(6,13) 

17)=(R*SIGINV(3))/(4.0*W(5)) 

15)=-U6/W(5)+DT3(6,17) 

16)=-U43/(2.0*W(5)> 

13)=(R*SIGINV(5))/(8.0*U(5)) 

14)=-DT3(7,13) 

15)=-U7/W(5) 

20)=-U73/(2.0*W(5)) 

21)=DT3(4,21) 

22)=<R*SIGINV(4))/(2.0*W(5)) 

11)=DT3(4,21) 

12)=-DT3(4f21) 

13)=DT3(7,22)/2.0 

14)=-DT3(8,13) 

15)=-U8/W(5)+DT3(7,13)*2.0 

16)=DT3(7,20) 

17)=DT3(7,13)*2.0 

11)=DT3(4,11) 

12)=-DT3(9,11) 
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DT3(9,13)=DT3(7,13) 

DT3(9,14)=-DT3(9,13) 

DT3(9,17)=(R*SIGINV(6)>/(8.0*U(5)) 

DT3(9,15)=-U9/W(5)+DT3(9,17) 

DT3(9,16)=-U93/(4.0*W(5)) 

RETURN 

END 


SUBROUTINE  MAKES 
CREATES  THE  MATRICES  S  AND  ST 


INCLUDE  'PMUCOM  DEF' 

INTEGER  I, J 

DO  10  1=1,9 
DO  20  J=1,6 
S(I,J)=0.0 
20     CONTINUE 
10  CONTINUE 

S(1,1)=1.0 

S(2,2)=0.5 

S<3,3)=0.5 

S<4,2)=0.5 

S(5,4>=1.0 

S(6,5)=0.5 

S(7,3)=0.5 

S(8,5)=0.5 

S<9,6)=1.0 

DO  30  1=1,6 
DO  40  J=1,9 

ST(I,J)=S(J,I) 
40     CONTINUE 
30  CONTINUE 

RETURN 

END 


SUBROUTINE  PMUS(INP,OUTP1 ,OUTP2) 

COMPUTES  THE  VALUES  OF  PMU  AND  PSIG 

INCLUDE  'PMUCOM  DEF' 
INCLUDE  'COM  DEF' 

REAL  INP(3,3),OUTP1(3,3),PRCO(9,9),TEMP(9,6),OUTP2(6,6) 
INTEGER  I, J 

CALL  TENSR(INP,INP,PROD,3,3,3,3) 
CALL  MXTMLT(PR0D,S,TEMP,9,9,6) 
CALL  MXTMLT(ST,TEMP,0UTP2,6,9,6) 
DO  10  1=1,6 
DO  20  J=1,6 

OUTP2(I,J)=(2.0*OUTP2(I,J))/(REAL(NUMTR)*NUMSIM) 
20     CONTINUE 
10  CONTINUE 

DO  30  1=1,3 
DO  40  J=1,3 

OUTP1(I,J)=INP(I,J)/REAL(NUMTR) 
40     CONTINUE 
30  CONTINUE 
RETURN 
END 
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SUBROUTINE  TENSR(A,B,C,A1 ,A2,B1 ,B2) 

COMPUTES  THE  TENSOR  PRODUCT  OF  A  AND  B  OUPTUT  TO  C 

INTEGER  I,J,K,L,A1,A2,B1,B2 
REAL  A(A1,A2),B(B1,B2),C(A1*B1,A2*B2) 
DO  10  1=1, A1 
DO  20  J=1,A2 
DO  30  K=1,B1 
DO  40  L=1,B2 

C(K+(I-1)*B1,L+(J-1)*B2)=A(I,J)*B(K,L) 
40  CONTINUE 

30        CONTINUE 
20     CONTINUE 
10  CONTINUE 
RETURN 
END 


SUBROUTINE  BNDEST(SEPM,PMU,PSIG,SEPMU,SEPSIG,CI ) 

CALCULATES  THE  FIRST  AND  SECOND  ORDER  VARIANCE  ESTIMATES  AND 
CORRESPONDING  CONFIDENCE  INTERVAL  SIZES 

INTEGER  I, J 

REAL  TEMP1(3,3),TEMP2(6,6),TEMP3(3f6),TEMP4(6,3),SSQMU(9), 
C     SSQSIG(36),SSQMSG(18),SSQSGM(18),SEPM(9,9),PSIGSQ(36,36), 
C     PSIG(6,6),SEPMU(3),SEPSIG(6),TEMP(36),PMUSQ(9,9),PMU(3,3), 
C     PMUSIG<18,18),PSIGMU(18,18),CI(2),STPMU(9),STPSIG(36) 
C     ,BND,BND1,BND2,BND3 
DO  10  1=1,3 
DO  20  J=1,3 

TEMP1(I,J)=SEPM(I,J) 
20     CONTINUE 

DO  30  J=1,6 

TEMP3(I,J)=SEPM(I,J+3) 
30     CONTINUE 
10  CONTINUE 

CALL  STRING(TEMP1,SSQMU,3,3) 
CALL  STRING(TEMP3,SSQMSG,3,6) 
DO  40  1=1,6 
DO  50  J=1,3 

TEMP4(I,J)=SEPM(I+3,J) 
50     CONTINUE 

DO  60  J=1,6 

TEMP2(I,J)=SEPM(I+3,J+3) 
60     CONTINUE 
40  CONTINUE 

CALL  STRING(TEMP4,SSQSGM,6,3) 
CALL  STRING(TEMP2,SSQSIG,6,6) 
CALL  MXTMLT(PMU,SEPMU,TEMP,3,3,1) 
CALL  MXTMLT(SEPMU,TEMP,BND, 1 ,3, 1 ) 
CALL  MXTMLT(PSIG,SEPSIG,TEMP,6,6,1) 
CALL  MXTMLT(SEPSIG,TEMP,BND1 , 1 ,6, 1 ) 
BND1=BND1+BND 

CALL  TENSR(PMU,PMU,PMUSQ,3,3,3,3) 
CALL  MXTMLT(PMUSQ,SSQMU,TEMP,9,9,1) 
CALL  MXTMLT(SSOMU, TEMP, BND, 1,9,1) 
CALL  TENSR(PSIG,PSIG,PSIGSQ,6,6,6,6) 
CALL  MXTMLT(PSIGSQ,SSQSIG, TEMP, 36, 36,1) 
CALL  MXTMLT(SSQSIG, TEMP, BND2, 1,36,1) 
BND2=(BND+BND2)/2.0 
CALL  TENSR(PMU,PSIG,PMUSIG,3,3,6,6) 
CALL  MXTMLT(PMUSIG,SSQMSG,TEMP,18,18,1) 
CALL  MXTMLT(SSQMSG, TEMP, BND, 1 ,18,1) 
CALL  TENSR(PSIG,PMU,PSIGMU,6,6,3,3) 
CALL  MXTMLT(PSIGMU,SSQSGM, TEMP, 18, 18,1) 
CALL  MXTMLT(SSQSGM,TEMP,BND3, 1,18,1) 
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BND2=BND2+((BND+BND3)/2.0) 

CALL  STRING(PMU,STPMU,3,3) 

CALL  STRING(PSIG,STPSIG,6,6) 

CALL  MXTMLT(SSQMU,STPMU,BND, 1,9,1) 

CALL  MXTMLT(SSQSIG,STPSIG,BND3, 1,36,1) 

BND2=BND2+( ( (BND+BND3)**2)/4 .0) 

BND=BND1+BND2 

CI(1)=1.96*S0RT(BND1) 

CI(2)=1.96*SQRT(BND) 

RETURN 

END 


SUBROUTINE  STRING(A1,A2,N1 ,N2) 

STRINGS  OUT  THE  MATRIX  A1  AS  A  VECTOR  A2 

REAL  A1(N1,N2),A2(N1*N2) 
INTEGER  I,N1,N2,J 
DO  10  1=1, N1 
DO  20  J=1,N2 

A2((I-1)*N2+J)=A1(I,J) 
20     CONTINUE 
10  CONTINUE 
RETURN 
END 


REAL  FUNCTION  HX(X,Y) 

REAL  X,Y,SPHER(3),TEMP 

INCLUDE  'COM  DEF' 

SPHER(1)=R*SIN(X)*COS(Y)-MU(1) 

SPHER(2)=R*SIN(X)*SIN(Y)-MU(2) 

SPHER(3)=R*COS(X)-MU(3) 

TEMP=-0.5*(((SPHER(1)**2)*SIGINV(1))+(2.0*SPHER(1)*SPHER(2) 
C  *SIGINV(2))+(2.0*SPHER(1)*SPHER(3)*SIGINV(3)) 

C  +<<SPHER(2)**2)*SIGINV(4))+(2.0*SPHER(2)*SPHER(3) 

C  *SIGINV(5))+((SPHER(3)**2)*SIGINV(6))) 

IF  (TEMP  .GT.  -27.0  )  THEN 
HX=EXP(TEMP)/DENOM 

ELSE 

HX=0.0 

ENDIF 

RETURN 

END 

REAL  FUNCTION  CC1111(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

CC1111=COS(X)*COS(Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CC1311(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

CC1311=COS(3.0*X)*COS(Y)*HX(X,Y) 

RETURN 

END 


REAL  FUNCTION  CC3111(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

CC3111=(COS(X)**3)*COS(Y)*HX(X,Y) 

RETURN 

END 
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REAL  FUNCTION  CS1111(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

CS1111=COS(X)*SIN(Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CS1311(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

CS1311=COS(3.0*X)*SIN(Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CS3111(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

CS3111=(COS(X)**3)*SIN(Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC1100(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

SC1100=SIN(X)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC1200(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

SC1200=SIN(2.0*X)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC1300(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

SC1300=SIN(3.0*X)*HX(X,Y) 

RETURN 

END 


REAL  FUNCTION  SC2111(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

SC2111=(SIN(X)**2)*COS(Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC3121(X,Y> 

EXTERNAL  HX 

REAL  X,Y 

SC3121=(SIN(X)**3)*(COS(Y)**2)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SS2111(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

SS2111=(SIN(X)**2)*SIN(Y)*HX(X,Y) 

RETURN 

END 
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REAL  FUNCTION  SS3112(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

SS3112=(SIN(X)**3)*SIN(2.0*Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SS3121(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

SS3121=(SIN(X)**3)*(SIN(Y)**2)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CC1113(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

CC1113=COS(X)*COS(3.0*Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CC1131(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

CC1131=C0S(X)*(C0S(Y)**3)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CC1211(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

CC1211=COS(2.0*X)*COS(Y)*HX(X,Y> 

RETURN 

END 

REAL  FUNCTION  CC1331(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

CC1331=COS(3.0*X)*(COS(Y)**3)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CC1411(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

CCH11=COS(4.0*X)*COS(Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CC1511(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

CC1511=COS(5.0*X)*COS(Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CC1531(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

CC1531=COS(5.0*X)*(COS(Y)**3)*HX(X,Y) 

RETURN 

END 
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REAL  FUNCTION  CC2111(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

CC2111=(COS(X)**2)*COS(Y)*HX<X,Y) 

RETURN 

END 

REAL  FUNCTION  CC3113(X,Y> 

EXTERNAL  HX 

REAL  X,Y 

CC3113=(COS(X)**3)*COS(3.0*Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CC3131(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

CC3131=(COS(X)**3)*(COS(Y)**3)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CC5111(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

CC5111=(COS(X)**5)*COS(Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CC5113(XfY) 

EXTERNAL  HX 

REAL  X,Y 

CC5113=(COS(X)**5)*COS(3.0*Y)*HX(X#Y) 

RETURN 

END 

REAL  FUNCTION  CC5131(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

CC5131=(COS(X)**5)*(COS(Y)**3)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CS1113(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

CS1113=COS(X)*SIN(3.0*Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CS1131(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

CS1131=COS(X)*(SIN(Y>**3)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CS1211(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

CS1211=COS(2.0*X)*SIN(Y)*HX(X,Y) 

RETURN 

END 
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REAL  FUNCTION  CS1331(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

CS1331=COS(3.0*X)*(SIN(Y)**3)*HX(X,Y> 

RETURN 

END 

REAL  FUNCTION  CSU11(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

CSH11=COS(4.0*X)*SIN(Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CS1511(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

CS1511=COS(5.0*X)*SIN(Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CS1531(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

CS1531=COS(5.0*X)*(SIN(Y)**3)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CS2111(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

CS2111=(COS(X)**2)*SIN(Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CS3113(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

CS3113=(COS(X)**3)*SIN(3.0*Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CS3131(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

CS3131=(COS(X)**3)*(SIN(Y)**3)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CS5111(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

CS5111=(COS(X)**5)*SIN(Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CS5113(X,Y) 

EXTERNAL  HX 

REAL  XfY 

CS5113=(COS(X)**5)*SIN(3.0*Y)*HX(X,Y) 

RETURN 

END 
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REAL  FUNCTION  CS5131(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

CS5131=(COS(X)**5)*(SIN(Y)**3)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC1121(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

SC1121=SIN(X)*(C0S(Y)**2)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC1221(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

SC1221=SIN(2.0*X)*(COS(Y)**2)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC1321(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

SC1321=SIN(3.0*X)*(COS(Y)**2)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC1400(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

SCH00=SIN(4.0*X)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC1421(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

SCH21=SIN(A.0*X)*(COS(Y)**2)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC1500(X,Y> 

EXTERNAL  HX 

REAL  X,Y 

SC1500=SIN(5.0*X)*HX(X,Y) 

RETURN 

END 

REAL    FUNCTION   SC152KX.Y) 

EXTERNAL  HX 

REAL  X,Y 

SC1521=SIN(5.0*X)*(COS(Y)**2)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC3100(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

SC3100=(SIN(X)**3)*HX(X,Y) 

RETURN 

END 
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REAL  FUNCTION  SC4111(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

SC4111=(SIN(X)**4)*COS(Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC4113(X,Y> 

EXTERNAL  HX 

REAL  X,Y 

SC4113=(SIN(X)**4)*COS(3.0*Y)*HX(X,Y> 

RETURN 

END 

REAL  FUNCTION  SC4131(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

SC4131=(SIN(X)**4)*<COS(Y)**3>*HX<X,Y) 

RETURN 

END 

REAL  FUNCTION  SC5112(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

SC5112=(SIN(X)**5)*COS(2.0*Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC51K(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

SC5114=(SIN(X)**5)*COS(4.0*Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC5121(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

SC5121=(SIN(X)**5)*(COS(Y)**2)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC5141(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

SC5141=(SIN(X)**5)*(COS(Y)**4)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC1212(X,Y) 

EXTERNAL  HX 

REAL  XfY 

SC1212=SIN(2.0*X)*COS(2.0*Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC1412(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

SC1412=SIN(4.0*X)*COS(2.0*Y)*HX(X,Y) 

RETURN 

END 
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REAL  FUNCTION  SS1112(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

SS1112=SIN(X)*SIN(2.0*Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SS1121(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

SS1121=SIN(X)*(SIN(Y)**2)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SS1212(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

SS1212=SIN(2.0*X)*SIN(2.0*Y)*HX(X,Y) 

RETURN 

END 

REAL   FUNCTION   SS122KX.Y) 

EXTERNAL  HX 

REAL  X,Y 

SS1221=SIN(2.0*X)*(SIN(Y)**2)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SS1312(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

SS1312=SIN(3.0*X)*SIN(2.0*Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SS1321(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

SS1321=SIN(3.0*X)*(SIN(Y)**2)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SSK12(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

SSK12=SIN(4.0*X)*SIN(2.0*Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SS1421(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

SS1421=SIN(4.0*X)*(SIN(Y)**2)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SS1512(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

SS1512=SIN(5.0*X)*SIN(2.0*Y)*HX(X,Y) 

RETURN 

END 


87 


REAL  FUNCTION  SS1521(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

SS1521=SIN(5.0*X)*(SIN(Y)**2)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SS411KX.Y) 

EXTERNAL  HX 

REAL  X,Y 

SS4111=(SIN(X)**4)*SIN(Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SS4113(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

SS4113=(SIN(X)**4)*SIN<3.0*Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SS413KX.Y) 

EXTERNAL  HX 

REAL  X,Y 

SS4131=(SIN(X)**4)*(SIN(Y)**3)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SS5112(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

SS5112=(SIN(X)**5)*SIN(2.0*Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SS5114(X,Y) 

EXTERNAL  HX 

REAL  X,Y 

SS5114=(SIN(X)**5)*SIN(4.0*Y)*HX(X,Y) 

RETURN 

END 

REAL   FUNCTION  SS512KX.Y) 

EXTERNAL  HX 

REAL  X,Y 

SS5121=(SIN(X)**5)*(SIN(Y)**2)*HX(X,Y) 

RETURN 

END 

REAL   FUNCTION   SS514KX.Y) 

EXTERNAL   HX 

REAL  X,Y 

SS5141=(SIN(X)**5)*(SIN(Y)**4)*HX(X,Y) 

RETURN 

END 


COM  DEF 

REAL  R,SIG(6),SIGINV(6),DETfDENOM,MU(3),PI,A,B,ERRABS,ERRREL 
C  ,MSIG(3/3),MUSET(3),SIGSET(3,3),SEPSET 

INTEGER   IRULE.ISEED.NUMTR 

COMMON/COMI/R.SIG.SIGINV.DET.DENOM.MU.ERRABS.ERRREL.MSIG, 
C  MUSET.SIGSET.SEPSET.IRULE.NUMTR.ISEED 

PARAMETER(PI=3. 1415927) 
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WVEC  DEF 

REAL  W(12),CC(17),CS(18),SC(19),SS(18) 

COMMON/UCC/W,CC,CS,SC,SS 


PMUCOM  DEF 

REAL  S(9,6),ST(6,9),NUMS1M 

COMMON/COM3/S,ST,NUMSIM 


SEPCOM  DEF 

REAL  CXX,CYY,CZZ,SQPI INV,QUPITW,CSANT,C, RADIUS, MN(3) 
COMMON/COM2/CXX , CY Y , CZZ , CSANT , C , RAD  I US , MN 
PARAMETER(SQPIINV=0.56A1896,QUPITW=0. 05626977) 
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